| 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 // $Id: G4PropagatorInField.cc,v 1.42 2008/01/24 08:54:01 gcosmo Exp $
>> 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
>> 29 //
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> <<
39 <<
40 #include "G4PropagatorInField.hh" 42 #include "G4PropagatorInField.hh"
41 #include "G4ios.hh" 43 #include "G4ios.hh"
42 #include "G4SystemOfUnits.hh" << 44 #include <iomanip>
>> 45
43 #include "G4ThreeVector.hh" 46 #include "G4ThreeVector.hh"
44 #include "G4Material.hh" <<
45 #include "G4VPhysicalVolume.hh" 47 #include "G4VPhysicalVolume.hh"
46 #include "G4Navigator.hh" 48 #include "G4Navigator.hh"
47 #include "G4GeometryTolerance.hh" 49 #include "G4GeometryTolerance.hh"
48 #include "G4VCurvedTrajectoryFilter.hh" 50 #include "G4VCurvedTrajectoryFilter.hh"
49 #include "G4ChordFinder.hh" 51 #include "G4ChordFinder.hh"
50 #include "G4MultiLevelLocator.hh" <<
51 52
52 << 53 ///////////////////////////////////////////////////////////////////////////
53 // ------------------------------------------- <<
54 // Constructors and destructor <<
55 // 54 //
56 G4PropagatorInField::G4PropagatorInField( G4Na << 55 // Constructors and destructor
57 G4Fi << 56
58 G4VI << 57 G4PropagatorInField::G4PropagatorInField( G4Navigator *theNavigator,
>> 58 G4FieldManager *detectorFieldMgr )
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 fParticleIsLooping(false),
>> 65 fVerboseLevel(0),
>> 66 fMax_loop_count(1000),
>> 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 )
64 { 72 {
65 fEpsilonStep = (fDetectorFieldMgr != nullptr << 73 if(fDetectorFieldMgr) { fEpsilonStep = fDetectorFieldMgr->GetMaximumEpsilonStep();}
66 ? fDetectorFieldMgr->GetMaximum << 74 else { fEpsilonStep= 1.0e-5; }
67 << 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 85
73 fLargestAcceptableStep = 100.0 * meter; // << 86 // In case of too slow progress in finding Intersection Point
74 fMaxStepSizeMultiplier= 0.1 ; // 0.1 in << 87 // intermediates Points on the Track must be stored.
75 fMinBigDistance= 100. * CLHEP::mm; << 88 // Initialise the array of Pointers [max_depth+1] to do this
76 #ifdef G4DEBUG_FIELD << 89
77 G4cout << " PiF: Zero Step Threshold set to << 90 G4ThreeVector zeroV(0.0,0.0,0.0);
78 << fZeroStepThreshold / millimeter << 91 for (G4int idepth=0; idepth<max_depth+1; idepth++ )
79 << " mm." << G4endl; <<
80 G4cout << " PiF: Value of kCarTolerance = <<
81 << kCarTolerance / millimeter <<
82 << " mm. " << G4endl; <<
83 fVerboseLevel = 2; <<
84 fVerbTracePiF = true; <<
85 #endif <<
86 <<
87 // Defining Intersection Locator and his par <<
88 if ( vLocator == nullptr ) <<
89 { <<
90 fIntersectionLocator = new G4MultiLevelLoc <<
91 fAllocatedLocator = true; <<
92 } <<
93 else <<
94 { 92 {
95 fIntersectionLocator = vLocator; << 93 ptrInterMedFT[ idepth ] = new G4FieldTrack( zeroV, zeroV, 0., 0., 0., 0.);
96 fAllocatedLocator = false; <<
97 } 94 }
98 RefreshIntersectionLocator(); // Copy all <<
99 } 95 }
100 96
101 // ------------------------------------------- <<
102 // <<
103 G4PropagatorInField::~G4PropagatorInField() 97 G4PropagatorInField::~G4PropagatorInField()
104 { 98 {
105 if(fAllocatedLocator) { delete fIntersecti << 99 for ( G4int idepth=0; idepth<max_depth+1; idepth++)
>> 100 {
>> 101 delete ptrInterMedFT[idepth];
>> 102 }
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 { <<
129 GetChordFinder()->OnComputeStep(&pFieldTrack <<
130 const G4double deltaChord = GetChordFinder() <<
131 <<
132 // If CurrentProposedStepLength is too small 116 // If CurrentProposedStepLength is too small for finding Chords
133 // then return with no action (for now - TOD 117 // then return with no action (for now - TODO: some action)
134 // 118 //
135 const char* methodName = "G4PropagatorInFiel << 119 if(CurrentProposedStepLength<kCarTolerance)
136 if (CurrentProposedStepLength<kCarTolerance) <<
137 { 120 {
138 return kInfinity; 121 return kInfinity;
139 } 122 }
140 123
141 // Introducing smooth trajectory display (ja 124 // Introducing smooth trajectory display (jacek 01/11/2002)
142 // 125 //
143 if (fpTrajectoryFilter != nullptr) << 126 if (fpTrajectoryFilter)
144 { 127 {
145 fpTrajectoryFilter->CreateNewTrajectorySeg 128 fpTrajectoryFilter->CreateNewTrajectorySegment();
146 } 129 }
147 130
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 131 // Parameters for adaptive Runge-Kutta integration
170 132
171 G4double h_TrialStepSize; // 1st Step << 133 G4double h_TrialStepSize; // 1st Step Size
172 G4double TruePathLength = CurrentProposedSte << 134 G4double TruePathLength = CurrentProposedStepLength;
173 G4double StepTaken = 0.0; << 135 G4double StepTaken = 0.0;
174 G4double s_length_taken, epsilon; << 136 G4double s_length_taken, epsilon ;
175 G4bool intersects; << 137 G4bool intersects;
176 G4bool first_substep = true; << 138 G4bool first_substep = true;
177 139
178 G4double NewSafety; << 140 G4double NewSafety;
179 fParticleIsLooping = false; 141 fParticleIsLooping = false;
180 142
181 // If not yet done, 143 // If not yet done,
182 // Set the field manager to the local one 144 // Set the field manager to the local one if the volume has one,
183 // or to the global one 145 // or to the global one if not
184 // 146 //
185 if( !fSetFieldMgr ) << 147 if( !fSetFieldMgr ) fCurrentFieldMgr= FindAndSetFieldManager( pPhysVol );
186 { << 148 // For the next call, the field manager must again be set
187 fCurrentFieldMgr = FindAndSetFieldManager( << 149 fSetFieldMgr= false;
188 } <<
189 fSetFieldMgr = false; // For next call, the <<
190 150
191 G4FieldTrack CurrentState(pFieldTrack); << 151 GetChordFinder()->SetChargeMomentumMass(fCharge, fInitialMomentumModulus, fMass);
192 G4FieldTrack OriginalState = CurrentState; << 152
>> 153 G4FieldTrack CurrentState(pFieldTrack);
>> 154 G4FieldTrack OriginalState = CurrentState;
193 155
194 // If the Step length is "infinite", then an 156 // If the Step length is "infinite", then an approximate-maximum Step
195 // length (used to calculate the relative ac << 157 // length (used to calculate the relative accuracy) must be guessed.
196 // 158 //
197 if( CurrentProposedStepLength >= fLargestAcc 159 if( CurrentProposedStepLength >= fLargestAcceptableStep )
198 { 160 {
199 G4ThreeVector StartPointA, VelocityUnit; 161 G4ThreeVector StartPointA, VelocityUnit;
200 StartPointA = pFieldTrack.GetPosition(); 162 StartPointA = pFieldTrack.GetPosition();
201 VelocityUnit = pFieldTrack.GetMomentumDir( 163 VelocityUnit = pFieldTrack.GetMomentumDir();
202 164
203 G4double trialProposedStep = fMaxStepSizeM << 165 G4double trialProposedStep = 1.e2 * ( 10.0 * cm +
204 fNavigator->GetWorldVolume()->GetLogical 166 fNavigator->GetWorldVolume()->GetLogicalVolume()->
205 GetSolid()->DistanceToOut(St 167 GetSolid()->DistanceToOut(StartPointA, VelocityUnit) );
206 CurrentProposedStepLength = std::min( tria << 168 CurrentProposedStepLength= std::min( trialProposedStep,
207 fLar << 169 fLargestAcceptableStep );
208 } 170 }
209 epsilon = fCurrentFieldMgr->GetDeltaOneStep( << 171 epsilon = GetDeltaOneStep() / CurrentProposedStepLength;
>> 172 // G4double raw_epsilon= epsilon;
210 G4double epsilonMin= fCurrentFieldMgr->GetMi 173 G4double epsilonMin= fCurrentFieldMgr->GetMinimumEpsilonStep();
211 G4double epsilonMax= fCurrentFieldMgr->GetMa << 174 G4double epsilonMax= fCurrentFieldMgr->GetMaximumEpsilonStep();;
212 if( epsilon < epsilonMin ) { epsilon = epsi << 175 if( epsilon < epsilonMin ) epsilon = epsilonMin;
213 if( epsilon > epsilonMax ) { epsilon = epsi << 176 if( epsilon > epsilonMax ) epsilon = epsilonMax;
214 SetEpsilonStep( epsilon ); 177 SetEpsilonStep( epsilon );
215 178
216 // Values for Intersection Locator has to be << 179 // G4cout << "G4PiF: Epsilon of current step - raw= " << raw_epsilon
217 // case that CurrentFieldManager has changed << 180 // << " final= " << epsilon << G4endl;
218 // <<
219 RefreshIntersectionLocator(); <<
220 181
221 // Shorten the proposed step in case of earl << 182 // Shorten the proposed step in case of earlier problems (zero steps)
222 // 183 //
223 if( fNoZeroStep > fActionThreshold_NoZeroSte 184 if( fNoZeroStep > fActionThreshold_NoZeroSteps )
224 { 185 {
225 G4double stepTrial; 186 G4double stepTrial;
226 187
227 stepTrial = fFull_CurveLen_of_LastAttempt; << 188 stepTrial= fFull_CurveLen_of_LastAttempt;
228 if( (stepTrial <= 0.0) && (fLast_ProposedS << 189 if( (stepTrial <= 0.0) && (fLast_ProposedStepLength > 0.0) )
229 { << 190 stepTrial= fLast_ProposedStepLength;
230 stepTrial = fLast_ProposedStepLength; <<
231 } <<
232 191
233 G4double decreaseFactor = 0.9; // Unused d 192 G4double decreaseFactor = 0.9; // Unused default
234 if( (fNoZeroStep < fSevereActionThreshol 193 if( (fNoZeroStep < fSevereActionThreshold_NoZeroSteps)
235 && (stepTrial > 100.0*fZeroStepThreshol << 194 && (stepTrial > 1000.0*kCarTolerance) )
236 { 195 {
237 // Attempt quick convergence << 196 // Ensure quicker convergence
238 // 197 //
239 decreaseFactor= 0.25; << 198 decreaseFactor= 0.1;
240 } 199 }
241 else 200 else
242 { 201 {
243 // We are in significant difficulties, p 202 // We are in significant difficulties, probably at a boundary that
244 // is either geometrically sharp or betw 203 // is either geometrically sharp or between very different materials.
245 // Careful decreases to cope with tolera << 204 // Careful decreases to cope with tolerance are required.
246 // 205 //
247 if( stepTrial > 100.0*fZeroStepThreshold << 206 if( stepTrial > 1000.0*kCarTolerance )
248 decreaseFactor = 0.35; // Try decr << 207 decreaseFactor = 0.25; // Try slow decreases
249 } else if( stepTrial > 30.0*fZeroStepThr << 208 else if( stepTrial > 100.0*kCarTolerance )
250 decreaseFactor= 0.5; // Try yet << 209 decreaseFactor= 0.5; // Try slower decreases
251 } else if( stepTrial > 10.0*fZeroStepThr << 210 else if( stepTrial > 10.0*kCarTolerance )
252 decreaseFactor= 0.75; // Try even 211 decreaseFactor= 0.75; // Try even slower decreases
253 } else { << 212 else
254 decreaseFactor= 0.9; // Try very 213 decreaseFactor= 0.9; // Try very slow decreases
255 } <<
256 } 214 }
257 stepTrial *= decreaseFactor; 215 stepTrial *= decreaseFactor;
258 216
259 #ifdef G4DEBUG_FIELD 217 #ifdef G4DEBUG_FIELD
260 if( fVerboseLevel > 2 << 218 PrintStepLengthDiagnostic(CurrentProposedStepLength, decreaseFactor,
261 || (fNoZeroStep >= fSevereActionThreshol << 219 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 220 #endif
275 if( stepTrial == 0.0 ) // Change to mak << 221 if( stepTrial == 0.0 )
276 { 222 {
277 std::ostringstream message; << 223 G4cout << " G4PropagatorInField::ComputeStep "
278 message << "Particle abandoned due to l << 224 << " Particle abandoned due to lack of progress in field."
279 << G4endl << 225 << G4endl
280 << " Properties : " << pFieldT << 226 << " Properties : " << pFieldTrack << " "
281 << " Attempting a zero step = << 227 << G4endl;
282 << " while attempting to progr << 228 G4cerr << " G4PropagatorInField::ComputeStep "
283 << " trial steps. Will abandon << 229 << " ERROR : attempting a zero step= " << stepTrial << G4endl
284 G4Exception(methodName, "GeomNav1002", << 230 << " while attempting to progress after " << fNoZeroStep
285 fParticleIsLooping = true; << 231 << " trial steps. Will abandon step." << G4endl;
286 return 0; // = stepTrial; << 232 fParticleIsLooping= true;
>> 233 return 0; // = stepTrial;
287 } 234 }
288 if( stepTrial < CurrentProposedStepLength 235 if( stepTrial < CurrentProposedStepLength )
289 { <<
290 CurrentProposedStepLength = stepTrial; 236 CurrentProposedStepLength = stepTrial;
291 } <<
292 } 237 }
293 fLast_ProposedStepLength = CurrentProposedSt 238 fLast_ProposedStepLength = CurrentProposedStepLength;
294 239
295 G4int do_loop_count = 0; 240 G4int do_loop_count = 0;
296 do // Loop checking, 07.10.2016, JA << 241 do
297 { 242 {
298 G4FieldTrack SubStepStartState = CurrentSt 243 G4FieldTrack SubStepStartState = CurrentState;
299 G4ThreeVector SubStartPoint = CurrentState 244 G4ThreeVector SubStartPoint = CurrentState.GetPosition();
300 << 245
301 if(!first_substep) << 246 if( !first_substep) {
302 { <<
303 if( fVerboseLevel > 4 ) <<
304 { <<
305 G4cout << " PiF: Calling Nav/Locate Gl <<
306 << G4endl; <<
307 } <<
308 fNavigator->LocateGlobalPointWithinVolum 247 fNavigator->LocateGlobalPointWithinVolume( SubStartPoint );
309 } 248 }
310 249
311 // How far to attempt to move the particle 250 // How far to attempt to move the particle !
312 // 251 //
313 h_TrialStepSize = CurrentProposedStepLengt 252 h_TrialStepSize = CurrentProposedStepLength - StepTaken;
314 253
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 254 // Integrate as far as "chord miss" rule allows.
326 // 255 //
327 s_length_taken = GetChordFinder()->Advance 256 s_length_taken = GetChordFinder()->AdvanceChordLimited(
328 CurrentState, 257 CurrentState, // Position & velocity
329 h_TrialStepSize, 258 h_TrialStepSize,
330 fEpsilonStep, 259 fEpsilonStep,
331 fPreviousSftOrigi 260 fPreviousSftOrigin,
332 fPreviousSafety ) << 261 fPreviousSafety
333 // CurrentState is now updated with the << 262 );
>> 263 // CurrentState is now updated with the final position and velocity.
334 264
335 fFull_CurveLen_of_LastAttempt = s_length_t 265 fFull_CurveLen_of_LastAttempt = s_length_taken;
336 266
337 G4ThreeVector EndPointB = CurrentState.Get << 267 G4ThreeVector EndPointB = CurrentState.GetPosition();
338 G4ThreeVector InterSectionPointE; << 268 G4ThreeVector InterSectionPointE;
339 G4double LinearStepLength; << 269 G4double LinearStepLength;
340 270
341 // Intersect chord AB with geometry 271 // Intersect chord AB with geometry
342 // <<
343 intersects= IntersectChord( SubStartPoint, 272 intersects= IntersectChord( SubStartPoint, EndPointB,
344 NewSafety, Lin << 273 NewSafety, LinearStepLength,
345 InterSectionPo 274 InterSectionPointE );
346 // E <- Intersection Point of chord AB a << 275 // E <- Intersection Point of chord AB and either volume A's surface
347 // or a << 276 // or a daughter volume's surface ..
348 277
349 if( first_substep ) << 278 if( first_substep ) {
350 { <<
351 currentSafety = NewSafety; 279 currentSafety = NewSafety;
352 } // Updating safety in other steps is pot 280 } // Updating safety in other steps is potential future extention
353 281
354 if( intersects ) 282 if( intersects )
355 { 283 {
356 G4FieldTrack IntersectPointVelct_G(Curr 284 G4FieldTrack IntersectPointVelct_G(CurrentState); // FT-Def-Construct
357 285
358 // Find the intersection point of AB tr 286 // Find the intersection point of AB true path with the surface
359 // of vol(A), if it exists. Start wit 287 // of vol(A), if it exists. Start with point E as first "estimate".
360 G4bool recalculatedEndPt = false; << 288 G4bool recalculatedEndPt= false;
361 << 289 G4bool found_intersection =
362 G4bool found_intersection = fIntersecti << 290 LocateIntersectionPoint( SubStepStartState, CurrentState,
363 EstimateIntersectionPoint( SubStepSta << 291 InterSectionPointE, IntersectPointVelct_G,
364 InterSecti << 292 recalculatedEndPt);
365 recalculat << 293 //G4cout<<"In Locate"<<recalculatedEndPt<<" and V"<<IntersectPointVelct_G.GetPosition()<<G4endl;
366 fPreviousS << 294 intersects = intersects && found_intersection;
367 intersects = found_intersection; << 295 if( found_intersection ) {
368 if( found_intersection ) <<
369 { <<
370 End_PointAndTangent= IntersectPointV 296 End_PointAndTangent= IntersectPointVelct_G; // G is our EndPoint ...
371 StepTaken = TruePathLength = Interse 297 StepTaken = TruePathLength = IntersectPointVelct_G.GetCurveLength()
372 - Origina << 298 - OriginalState.GetCurveLength();
373 } << 299 } else {
374 else << 300 // intersects= false; // "Minor" chords do not intersect
375 { << 301 if( recalculatedEndPt ){
376 // Either "minor" chords do not inte << 302 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 } 303 }
401 } 304 }
402 } 305 }
403 if( !intersects ) 306 if( !intersects )
404 { 307 {
405 StepTaken += s_length_taken; << 308 StepTaken += s_length_taken;
406 << 309 // For smooth trajectory display (jacek 01/11/2002)
407 if (fpTrajectoryFilter != nullptr) // Fo << 310 if (fpTrajectoryFilter) {
408 { <<
409 fpTrajectoryFilter->TakeIntermediatePo 311 fpTrajectoryFilter->TakeIntermediatePoint(CurrentState.GetPosition());
410 } 312 }
411 } 313 }
412 first_substep = false; 314 first_substep = false;
413 315
414 #ifdef G4DEBUG_FIELD 316 #ifdef G4DEBUG_FIELD
415 if( fNoZeroStep > fActionThreshold_NoZeroS << 317 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 318 printStatus( SubStepStartState, // or OriginalState,
423 CurrentState, CurrentPropos << 319 CurrentState, CurrentProposedStepLength,
424 NewSafety, do_loop_count, p << 320 NewSafety, do_loop_count, pPhysVol );
425 } 321 }
426 if( (fVerboseLevel > 1) && (do_loop_count << 322 #endif
427 { << 323 #ifdef G4VERBOSE
428 if( do_loop_count == fMax_loop_count-9 ) << 324 if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 )) {
429 { << 325 if( do_loop_count == fMax_loop_count-9 ){
430 G4cout << " G4PropagatorInField::Compu << 326 G4cout << "G4PropagatorInField::ComputeStep "
431 << " Difficult track - taking << 327 << " Difficult track - taking many sub steps." << G4endl;
432 printStatus( SubStepStartState, SubSte <<
433 NewSafety, 0, pPhysVol ); <<
434 } 328 }
435 printStatus( SubStepStartState, CurrentS 329 printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength,
436 NewSafety, do_loop_count, p 330 NewSafety, do_loop_count, pPhysVol );
437 } 331 }
438 #endif 332 #endif
439 333
440 ++do_loop_count; << 334 do_loop_count++;
441 335
442 } while( (!intersects ) 336 } while( (!intersects )
443 && (!fParticleIsLooping) <<
444 && (StepTaken + kCarTolerance < Curren 337 && (StepTaken + kCarTolerance < CurrentProposedStepLength)
445 && ( do_loop_count < fMax_loop_count ) 338 && ( do_loop_count < fMax_loop_count ) );
446 339
447 if( do_loop_count >= fMax_loop_count << 340 if( do_loop_count >= fMax_loop_count )
448 && (StepTaken + kCarTolerance < CurrentPro <<
449 { 341 {
450 fParticleIsLooping = true; 342 fParticleIsLooping = true;
>> 343
>> 344 if ( fVerboseLevel > 0 ){
>> 345 G4cout << "G4PropagateInField: Killing looping particle "
>> 346 // << " of " << energy << " energy "
>> 347 << " after " << do_loop_count << " field substeps "
>> 348 << " totaling " << StepTaken / mm << " mm " ;
>> 349 if( pPhysVol )
>> 350 G4cout << " in the volume " << pPhysVol->GetName() ;
>> 351 else
>> 352 G4cout << " in unknown or null volume. " ;
>> 353 G4cout << G4endl;
>> 354 }
451 } 355 }
452 if ( ( fParticleIsLooping ) && (fVerboseLeve << 356
453 { <<
454 ReportLoopingParticle( do_loop_count, Step <<
455 CurrentProposedStep <<
456 CurrentState.GetMom <<
457 } <<
458 <<
459 if( !intersects ) 357 if( !intersects )
460 { 358 {
461 // Chord AB or "minor chords" do not inter 359 // Chord AB or "minor chords" do not intersect
462 // B is the endpoint Step of the current S 360 // B is the endpoint Step of the current Step.
463 // 361 //
464 End_PointAndTangent = CurrentState; 362 End_PointAndTangent = CurrentState;
465 TruePathLength = StepTaken; // Original << 363 TruePathLength = StepTaken;
466 <<
467 // Tried the following to avoid potential <<
468 // - but has issues... Suppressing this ch <<
469 // TruePathLength = CurrentProposedStepLen <<
470 } 364 }
471 fLastStepInVolume = intersects; <<
472 365
473 // Set pFieldTrack to the return value 366 // Set pFieldTrack to the return value
474 // 367 //
475 pFieldTrack = End_PointAndTangent; 368 pFieldTrack = End_PointAndTangent;
476 369
477 #ifdef G4VERBOSE 370 #ifdef G4VERBOSE
478 // Check that "s" is correct 371 // Check that "s" is correct
479 // 372 //
480 if( std::fabs(OriginalState.GetCurveLength() 373 if( std::fabs(OriginalState.GetCurveLength() + TruePathLength
481 - End_PointAndTangent.GetCurveLength()) 374 - End_PointAndTangent.GetCurveLength()) > 3.e-4 * TruePathLength )
482 { 375 {
483 std::ostringstream message; << 376 G4cerr << " ERROR - G4PropagatorInField::ComputeStep():" << G4endl
484 message << "Curve length mis-match between << 377 << " Curve length mis-match, is advancement wrong ? " << G4endl;
485 << "and proposed endpoint of propa << 378 G4cerr << " The curve length of the endpoint should be: "
486 << " The curve length of the endp << 379 << OriginalState.GetCurveLength() + TruePathLength << G4endl
487 << OriginalState.GetCurveLength() << 380 << " and it is instead: "
488 << " and it is instead: " << 381 << End_PointAndTangent.GetCurveLength() << "." << G4endl
489 << End_PointAndTangent.GetCurveLen << 382 << " A difference of: "
490 << " A difference of: " << 383 << OriginalState.GetCurveLength() + TruePathLength
491 << OriginalState.GetCurveLength() << 384 - End_PointAndTangent.GetCurveLength() << G4endl;
492 - End_PointAndTangent.GetCurveL << 385 G4cerr << " Original state= " << OriginalState << G4endl
493 << " Original state = " << Origin << 386 << " Proposed state= " << End_PointAndTangent << G4endl;
494 << " Proposed state = " << End_Po << 387 G4Exception("G4PropagatorInField::ComputeStep()", "IncorrectProposedEndPoint",
495 G4Exception(methodName, "GeomNav0003", Fat << 388 FatalException,
>> 389 "Curve length mis-match between original state and proposed endpoint of propagation.");
496 } 390 }
497 #endif 391 #endif
498 392
499 if( TruePathLength+kCarTolerance >= CurrentP << 393 // In particular anomalous cases, we can get repeated zero steps
500 { << 394 // In order to correct this efficiently, we identify these cases
501 fNoZeroStep = 0; << 395 // and only take corrective action when they occur.
502 } << 396 //
>> 397 if( TruePathLength < 0.5*kCarTolerance )
>> 398 fNoZeroStep++;
503 else 399 else
504 { << 400 fNoZeroStep = 0;
505 // In particular anomalous cases, we can << 401
506 // We identify these cases and take corre << 402 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; 403 fParticleIsLooping = true;
520 ReportStuckParticle( fNoZeroStep, Current << 404 G4cout << " WARNING - G4PropagatorInField::ComputeStep():" << G4endl
521 fFull_CurveLen_of_La << 405 << " Zero progress for " << fNoZeroStep << " attempted steps."
>> 406 << G4endl;
>> 407 if ( fVerboseLevel > 2 )
>> 408 G4cout << " Particle that is stuck will be killed." << G4endl;
522 fNoZeroStep = 0; 409 fNoZeroStep = 0;
523 } 410 }
524 << 411 // G4cout << "G4PropagatorInField returns " << TruePathLength << G4endl;
525 GetChordFinder()->SetDeltaChord(deltaChord); <<
526 return TruePathLength; 412 return TruePathLength;
527 } 413 }
528 414
529 // ------------------------------------------- << 415 // --------------------------------------------------------------------------
530 // Dumps status of propagator << 416 // G4bool
>> 417 // G4PropagatorInField::LocateIntersectionPoint(
>> 418 // const G4FieldTrack& CurveStartPointVelocity, // A
>> 419 // const G4FieldTrack& CurveEndPointVelocity, // B
>> 420 // const G4ThreeVector& TrialPoint, // E
>> 421 // G4FieldTrack& IntersectedOrRecalculated // Output
>> 422 // G4bool& recalculated) // Out
>> 423 // --------------------------------------------------------------------------
>> 424 //
>> 425 // Function that returns the intersection of the true path with the surface
>> 426 // of the current volume (either the external one or the inner one with one
>> 427 // of the daughters
>> 428 //
>> 429 // A = Initial point
>> 430 // B = another point
>> 431 //
>> 432 // Both A and B are assumed to be on the true path.
>> 433 //
>> 434 // E is the first point of intersection of the chord AB with
>> 435 // a volume other than A (on the surface of A or of a daughter)
>> 436 //
>> 437 // Convention of Use :
>> 438 // i) If it returns "true", then IntersectionPointVelocity is set
>> 439 // to the approximate intersection point.
>> 440 // ii) If it returns "false", no intersection was found.
>> 441 // The validity of IntersectedOrRecalculated depends on 'recalculated'
>> 442 // a) if latter is false, then IntersectedOrRecalculated is invalid.
>> 443 // b) if latter is true, then IntersectedOrRecalculated is
>> 444 // the new endpoint, due to a re-integration.
>> 445 // --------------------------------------------------------------------------
>> 446
>> 447 G4bool
>> 448 G4PropagatorInField::LocateIntersectionPoint(
>> 449 const G4FieldTrack& CurveStartPointVelocity, // A
>> 450 const G4FieldTrack& CurveEndPointVelocity, // B
>> 451 const G4ThreeVector& TrialPoint, // E
>> 452 G4FieldTrack& IntersectedOrRecalculatedFT, // Out: point found
>> 453 G4bool& recalculatedEndPoint) // Out:
>> 454 {
>> 455 // Find Intersection Point ( A, B, E ) of true path AB - start at E.
>> 456
>> 457 G4bool found_approximate_intersection = false;
>> 458 G4bool there_is_no_intersection = false;
>> 459
>> 460 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
>> 461 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
>> 462 G4ThreeVector CurrentE_Point = TrialPoint;
>> 463 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct
>> 464 G4double NewSafety= -0.0;
>> 465
>> 466 G4bool final_section= true; // Shows whether current section is last
>> 467 // (i.e. B=full end)
>> 468 G4bool first_section=true;
>> 469 recalculatedEndPoint= false;
>> 470
>> 471 G4bool restoredFullEndpoint= false;
>> 472
>> 473 G4int substep_no = 0;
>> 474
>> 475 // Limits for substep number
>> 476 //
>> 477 const G4int max_substeps= 10000; // Test 120 (old value 100 )
>> 478 const G4int warn_substeps= 1000; // 100
>> 479
>> 480 // Statistics for substeps
>> 481 //
>> 482 static G4int max_no_seen= -1;
>> 483 static G4int trigger_substepno_print= warn_substeps - 20 ;
>> 484
>> 485 //--------------------------------------------------------------------------
>> 486 // Algoritm for the case if progress in founding intersection is too slow.
>> 487 // Process is defined too slow if after N=param_substeps advances on the
>> 488 // path, it will be only 'fraction_done' of the total length.
>> 489 // In this case the remaining length is divided in two half and
>> 490 // the loop is restarted for each half.
>> 491 // If progress is still too slow, the division in two halfs continue
>> 492 // until 'max_depth'.
>> 493 //--------------------------------------------------------------------------
>> 494
>> 495 const G4int param_substeps=10; // Test value for the maximum number
>> 496 // of substeps
>> 497 const G4double fraction_done=0.3;
>> 498
>> 499 G4bool Second_half=false; // First half or second half of divided step
>> 500
>> 501 // We need to know this for the 'final_section':
>> 502 // real 'final_section' or first half 'final_section'
>> 503 // In algorithm it is considered that the 'Second_half' is true
>> 504 // and it becomes false only if we are in the first-half of level
>> 505 // depthness or if we are in the first section
>> 506
>> 507 G4int depth=0; // Depth counts how many subdivisions of initial step made
>> 508
>> 509 #ifdef G4DEBUG_FIELD
>> 510 static G4double tolerance= 1.0e-8;
>> 511 G4ThreeVector StartPosition= CurveStartPointVelocity.GetPosition();
>> 512 if( (TrialPoint - StartPosition).mag() < tolerance * mm )
>> 513 {
>> 514 G4cerr << "WARNING - G4PropagatorInField::LocateIntersectionPoint()"
>> 515 << G4endl
>> 516 << " Intermediate F point is on top of starting point A."
>> 517 << G4endl;
>> 518 G4Exception("G4PropagatorInField::LocateIntersectionPoint()",
>> 519 "IntersectionPointIsAtStart", JustWarning,
>> 520 "Intersection point F is exactly at start point A." );
>> 521 }
>> 522 #endif
>> 523
>> 524 // Intermediates Points on the Track = Subdivided Points must be stored.
>> 525 // Give the initial values to 'InterMedFt'
>> 526 // Important is 'ptrInterMedFT[0]', it saves the 'EndCurvePoint'
>> 527 //
>> 528 *ptrInterMedFT[0] = CurveEndPointVelocity;
>> 529 for (G4int idepth=1; idepth<max_depth+1; idepth++ )
>> 530 {
>> 531 *ptrInterMedFT[idepth]=CurveStartPointVelocity;
>> 532 }
>> 533
>> 534 // 'SubStartPoint' is needed to calculate the length of the divided step
>> 535 //
>> 536 G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
>> 537
>> 538 do
>> 539 {
>> 540 G4int substep_no_p = 0;
>> 541 G4bool sub_final_section = false; // the same as final_section,
>> 542 // but for 'sub_section'
>> 543 do // REPEAT param
>> 544 {
>> 545 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
>> 546 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
>> 547
>> 548 // F = a point on true AB path close to point E
>> 549 // (the closest if possible)
>> 550 //
>> 551 ApproxIntersecPointV = GetChordFinder()
>> 552 ->ApproxCurvePointV( CurrentA_PointVelocity,
>> 553 CurrentB_PointVelocity,
>> 554 CurrentE_Point,
>> 555 fEpsilonStep );
>> 556 // The above method is the key & most intuitive part ...
>> 557
>> 558 #ifdef G4DEBUG_FIELD
>> 559 if( ApproxIntersecPointV.GetCurveLength() >
>> 560 CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
>> 561 {
>> 562 G4cerr << "ERROR - G4PropagatorInField::LocateIntersectionPoint()"
>> 563 << G4endl
>> 564 << " Intermediate F point is more advanced than"
>> 565 << " endpoint B." << G4endl;
>> 566 G4Exception("G4PropagatorInField::LocateIntersectionPoint()",
>> 567 "IntermediatePointConfusion", FatalException,
>> 568 "Intermediate F point is past end B point" );
>> 569 }
>> 570 #endif
>> 571
>> 572 G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
>> 573
>> 574 // First check whether EF is small - then F is a good approx. point
>> 575 // Calculate the length and direction of the chord AF
>> 576 //
>> 577 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
>> 578
>> 579 if ( ChordEF_Vector.mag2() <= sqr(GetDeltaIntersection()) )
>> 580 {
>> 581 found_approximate_intersection = true;
>> 582
>> 583 // Create the "point" return value
>> 584 //
>> 585 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
>> 586 IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
>> 587
>> 588 // Note: in order to return a point on the boundary,
>> 589 // we must return E. But it is F on the curve.
>> 590 // So we must "cheat": we are using the position at point E
>> 591 // and the velocity at point F !!!
>> 592 //
>> 593 // This must limit the length we can allow for displacement!
>> 594 }
>> 595 else // E is NOT close enough to the curve (ie point F)
>> 596 {
>> 597 // Check whether any volumes are encountered by the chord AF
>> 598 // ---------------------------------------------------------
>> 599 // First relocate to restore any Voxel etc information
>> 600 // in the Navigator before calling ComputeStep()
>> 601 //
>> 602 fNavigator->LocateGlobalPointWithinVolume( Point_A );
>> 603
>> 604 G4ThreeVector PointG; // Candidate intersection point
>> 605 G4double stepLengthAF;
>> 606 G4bool Intersects_AF = IntersectChord( Point_A, CurrentF_Point,
>> 607 NewSafety, stepLengthAF,
>> 608 PointG );
>> 609 if( Intersects_AF )
>> 610 {
>> 611 // G is our new Candidate for the intersection point.
>> 612 // It replaces "E" and we will repeat the test to see if
>> 613 // it is a good enough approximate point for us.
>> 614 // B <- F
>> 615 // E <- G
>> 616
>> 617 CurrentB_PointVelocity = ApproxIntersecPointV;
>> 618 CurrentE_Point = PointG;
>> 619
>> 620 // By moving point B, must take care if current
>> 621 // AF has no intersection to try current FB!!
>> 622 //
>> 623 final_section= false;
>> 624
>> 625 #ifdef G4VERBOSE
>> 626 if( fVerboseLevel > 3 )
>> 627 {
>> 628 G4cout << "G4PiF::LI> Investigating intermediate point"
>> 629 << " at s=" << ApproxIntersecPointV.GetCurveLength()
>> 630 << " on way to full s="
>> 631 << CurveEndPointVelocity.GetCurveLength() << G4endl;
>> 632 }
>> 633 #endif
>> 634 }
>> 635 else // not Intersects_AF
>> 636 {
>> 637 // In this case:
>> 638 // There is NO intersection of AF with a volume boundary.
>> 639 // We must continue the search in the segment FB!
>> 640 //
>> 641 fNavigator->LocateGlobalPointWithinVolume( CurrentF_Point );
>> 642
>> 643 G4double stepLengthFB;
>> 644 G4ThreeVector PointH;
>> 645
>> 646 // Check whether any volumes are encountered by the chord FB
>> 647 // ---------------------------------------------------------
>> 648
>> 649 G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B,
>> 650 NewSafety, stepLengthFB,
>> 651 PointH );
>> 652 if( Intersects_FB )
>> 653 {
>> 654 // There is an intersection of FB with a volume boundary
>> 655 // H <- First Intersection of Chord FB
>> 656
>> 657 // H is our new Candidate for the intersection point.
>> 658 // It replaces "E" and we will repeat the test to see if
>> 659 // it is a good enough approximate point for us.
>> 660
>> 661 // Note that F must be in volume volA (the same as A)
>> 662 // (otherwise AF would meet a volume boundary!)
>> 663 // A <- F
>> 664 // E <- H
>> 665
>> 666 CurrentA_PointVelocity = ApproxIntersecPointV;
>> 667 CurrentE_Point = PointH;
>> 668 }
>> 669 else // not Intersects_FB
>> 670 {
>> 671 // There is NO intersection of FB with a volume boundary
>> 672
>> 673 if( final_section )
>> 674 {
>> 675 // If B is the original endpoint, this means that whatever
>> 676 // volume(s) intersected the original chord, none touch the
>> 677 // smaller chords we have used.
>> 678 // The value of 'IntersectedOrRecalculatedFT' returned is
>> 679 // likely not valid
>> 680
>> 681 // Check on real final_section or SubEndSection
>> 682 //
>> 683 if( ((Second_half)&&(depth==0)) || (first_section) )
>> 684 {
>> 685 there_is_no_intersection = true; // real final_section
>> 686 }
>> 687 else
>> 688 {
>> 689 // end of subsection, not real final section
>> 690 // exit from the and go to the depth-1 level
>> 691
>> 692 substep_no_p = param_substeps+2; // exit from the loop
>> 693
>> 694 // but 'Second_half' is still true because we need to find
>> 695 // the 'CurrentE_point' for the next loop
>> 696 //
>> 697 Second_half = true;
>> 698 sub_final_section = true;
>> 699
>> 700 }
>> 701 }
>> 702 else
>> 703 {
>> 704 // We must restore the original endpoint
>> 705
>> 706 CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
>> 707 CurrentB_PointVelocity = CurveEndPointVelocity;
>> 708 restoredFullEndpoint = true;
>> 709 }
>> 710 } // Endif (Intersects_FB)
>> 711 } // Endif (Intersects_AF)
>> 712
>> 713 // Ensure that the new endpoints are not further apart in space
>> 714 // than on the curve due to different errors in the integration
>> 715 //
>> 716 G4double linDistSq, curveDist;
>> 717 linDistSq = ( CurrentB_PointVelocity.GetPosition()
>> 718 - CurrentA_PointVelocity.GetPosition() ).mag2();
>> 719 curveDist = CurrentB_PointVelocity.GetCurveLength()
>> 720 - CurrentA_PointVelocity.GetCurveLength();
>> 721
>> 722 // Change this condition for very strict parameters of propagation
>> 723 //
>> 724 if( curveDist*curveDist*(1+2* fEpsilonStep ) < linDistSq )
>> 725 {
>> 726 // Re-integrate to obtain a new B
>> 727 //
>> 728 G4FieldTrack newEndPointFT=
>> 729 ReEstimateEndpoint( CurrentA_PointVelocity,
>> 730 CurrentB_PointVelocity,
>> 731 linDistSq, // to avoid recalculation
>> 732 curveDist );
>> 733 G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
>> 734 CurrentB_PointVelocity = newEndPointFT;
>> 735
>> 736 if( (final_section)&&(Second_half)&&(depth==0) ) // real final section
>> 737 {
>> 738 recalculatedEndPoint = true;
>> 739 IntersectedOrRecalculatedFT = newEndPointFT;
>> 740 // So that we can return it, if it is the endpoint!
>> 741 }
>> 742 }
>> 743 if( curveDist < 0.0 )
>> 744 {
>> 745 G4cerr << "ERROR - G4PropagatorInField::LocateIntersectionPoint()"
>> 746 << G4endl
>> 747 << " Error in advancing propagation." << G4endl;
>> 748 fVerboseLevel = 5; // Print out a maximum of information
>> 749 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
>> 750 -1.0, NewSafety, substep_no, 0 );
>> 751 G4cerr << " Point A (start) is " << CurrentA_PointVelocity
>> 752 << G4endl;
>> 753 G4cerr << " Point B (end) is " << CurrentB_PointVelocity
>> 754 << G4endl;
>> 755 G4cerr << " Curve distance is " << curveDist << G4endl;
>> 756 G4cerr << G4endl
>> 757 << "The final curve point is not further along"
>> 758 << " than the original!" << G4endl;
>> 759
>> 760 if( recalculatedEndPoint )
>> 761 {
>> 762 G4cerr << "Recalculation of EndPoint was called with fEpsStep= "
>> 763 << fEpsilonStep << G4endl;
>> 764 }
>> 765 G4cerr.precision(20);
>> 766 G4cerr << " Point A (Curve start) is " << CurveStartPointVelocity
>> 767 << G4endl;
>> 768 G4cerr << " Point B (Curve end) is " << CurveEndPointVelocity
>> 769 << G4endl;
>> 770 G4cerr << " Point A (Current start) is " << CurrentA_PointVelocity
>> 771 << G4endl;
>> 772 G4cerr << " Point B (Current end) is " << CurrentB_PointVelocity
>> 773 << G4endl;
>> 774 G4cerr << " Point S (Sub start) is " << SubStart_PointVelocity
>> 775 << G4endl;
>> 776 G4cerr << " Point E (Trial Point) is " << CurrentE_Point
>> 777 << G4endl;
>> 778 G4cerr << " Point F (Intersection) is " << ApproxIntersecPointV
>> 779 << G4endl;
>> 780 G4cerr << " LocateIntersection parameters are : Substep no= "
>> 781 << substep_no << G4endl;
>> 782 G4cerr << " Substep depth no= "<< substep_no_p << " Depth= "
>> 783 << depth << G4endl;
>> 784
>> 785 G4Exception("G4PropagatorInField::LocateIntersectionPoint()",
>> 786 "FatalError", FatalException,
>> 787 "Error in advancing propagation.");
>> 788 }
>> 789
>> 790 if(restoredFullEndpoint)
>> 791 {
>> 792 final_section = restoredFullEndpoint;
>> 793 restoredFullEndpoint = false;
>> 794 }
>> 795 } // EndIf ( E is close enough to the curve, ie point F. )
>> 796 // tests ChordAF_Vector.mag() <= maximum_lateral_displacement
>> 797
>> 798 #ifdef G4DEBUG_LOCATE_INTERSECTION
>> 799 if( substep_no >= trigger_substepno_print )
>> 800 {
>> 801 G4cout << "Difficulty in converging in "
>> 802 << "G4PropagatorInField::LocateIntersectionPoint():"
>> 803 << G4endl
>> 804 << " Substep no = " << substep_no << G4endl;
>> 805 if( substep_no == trigger_substepno_print )
>> 806 {
>> 807 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
>> 808 -1.0, NewSafety, 0, 0);
>> 809 }
>> 810 G4cout << " State of point A: ";
>> 811 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
>> 812 -1.0, NewSafety, substep_no-1, 0);
>> 813 G4cout << " State of point B: ";
>> 814 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
>> 815 -1.0, NewSafety, substep_no, 0);
>> 816 }
>> 817 #endif
>> 818
>> 819 substep_no++;
>> 820 substep_no_p++;
>> 821
>> 822 } while ( ( ! found_approximate_intersection )
>> 823 && ( ! there_is_no_intersection )
>> 824 && ( substep_no_p <= param_substeps) ); // UNTIL found or
>> 825 // failed param substep
>> 826 first_section = false;
>> 827
>> 828 if( (!found_approximate_intersection) && (!there_is_no_intersection) )
>> 829 {
>> 830 G4double did_len = std::abs( CurrentA_PointVelocity.GetCurveLength()
>> 831 - SubStart_PointVelocity.GetCurveLength());
>> 832 G4double all_len = std::abs( CurrentB_PointVelocity.GetCurveLength()
>> 833 - SubStart_PointVelocity.GetCurveLength());
>> 834
>> 835 G4double stepLengthAB;
>> 836 G4ThreeVector PointGe;
>> 837
>> 838 // Check if progress is too slow and if it possible to go deeper,
>> 839 // then halve the step if so
>> 840 //
>> 841 if( ( ( did_len )<fraction_done*all_len)
>> 842 && (depth<max_depth) && (!sub_final_section) )
>> 843 {
>> 844
>> 845 Second_half=false;
>> 846 depth++;
>> 847
>> 848 G4double Sub_len = (all_len-did_len)/(2.);
>> 849 G4FieldTrack start = CurrentA_PointVelocity;
>> 850 G4MagInt_Driver* integrDriver=GetChordFinder()->GetIntegrationDriver();
>> 851 integrDriver->AccurateAdvance(start, Sub_len, fEpsilonStep);
>> 852 *ptrInterMedFT[depth] = start;
>> 853 CurrentB_PointVelocity = *ptrInterMedFT[depth];
>> 854
>> 855 // Adjust 'SubStartPoint' to calculate the 'did_length' in next loop
>> 856 //
>> 857 SubStart_PointVelocity = CurrentA_PointVelocity;
>> 858
>> 859 // Find new trial intersection point needed at start of the loop
>> 860 //
>> 861 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
>> 862 G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();
>> 863
>> 864 fNavigator->LocateGlobalPointWithinVolume(Point_A);
>> 865 G4bool Intersects_AB = IntersectChord(Point_A, SubE_point,
>> 866 NewSafety, stepLengthAB, PointGe);
>> 867 if(Intersects_AB)
>> 868 {
>> 869 CurrentE_Point = PointGe;
>> 870 }
>> 871 else
>> 872 {
>> 873 // No intersection found for first part of curve
>> 874 // (CurrentA,InterMedPoint[depth]). Go to the second part
>> 875 //
>> 876 Second_half = true;
>> 877 }
>> 878 } // if did_len
>> 879
>> 880 if( (Second_half)&&(depth!=0) )
>> 881 {
>> 882 // Second part of curve (InterMed[depth],Intermed[depth-1]) )
>> 883 // On the depth-1 level normally we are on the 'second_half'
>> 884
>> 885 Second_half = true;
>> 886
>> 887 // Find new trial intersection point needed at start of the loop
>> 888 //
>> 889 SubStart_PointVelocity = *ptrInterMedFT[depth];
>> 890 CurrentA_PointVelocity = *ptrInterMedFT[depth];
>> 891 CurrentB_PointVelocity = *ptrInterMedFT[depth-1];
>> 892 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
>> 893 G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();
>> 894 fNavigator->LocateGlobalPointWithinVolume(Point_A);
>> 895 G4bool Intersects_AB = IntersectChord(Point_A, SubE_point, NewSafety,
>> 896 stepLengthAB, PointGe);
>> 897 if(Intersects_AB)
>> 898 {
>> 899 CurrentE_Point = PointGe;
>> 900 }
>> 901 else
>> 902 {
>> 903 final_section = true;
>> 904 }
>> 905 depth--;
>> 906 }
>> 907 } // if(!found_aproximate_intersection)
>> 908
>> 909 } while ( ( ! found_approximate_intersection )
>> 910 && ( ! there_is_no_intersection )
>> 911 && ( substep_no <= max_substeps) ); // UNTIL found or failed
>> 912
>> 913 if( substep_no > max_no_seen )
>> 914 {
>> 915 max_no_seen = substep_no;
>> 916 if( max_no_seen > warn_substeps )
>> 917 {
>> 918 trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps
>> 919 }
>> 920 }
>> 921
>> 922 if( ( substep_no >= max_substeps)
>> 923 && !there_is_no_intersection
>> 924 && !found_approximate_intersection )
>> 925 {
>> 926 G4cerr << "WARNING - G4PropagatorInField::LocateIntersectionPoint()"
>> 927 << G4endl
>> 928 << " Convergence is requiring too many substeps: "
>> 929 << substep_no << G4endl;
>> 930 G4cerr << " Abandoning effort to intersect. " << G4endl;
>> 931 G4cerr << " Information on start & current step follows in cout."
>> 932 << G4endl;
>> 933 G4cout << "WARNING - G4PropagatorInField::LocateIntersectionPoint()"
>> 934 << G4endl
>> 935 << " Convergence is requiring too many substeps: "
>> 936 << substep_no << G4endl;
>> 937 G4cout << " Found intersection = "
>> 938 << found_approximate_intersection << G4endl
>> 939 << " Intersection exists = "
>> 940 << !there_is_no_intersection << G4endl;
>> 941 G4cout << " Start and Endpoint of Requested Step:" << G4endl;
>> 942 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
>> 943 -1.0, NewSafety, 0, 0);
>> 944 G4cout << G4endl;
>> 945 G4cout << " 'Bracketing' starting and endpoint of current Sub-Step"
>> 946 << G4endl;
>> 947 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
>> 948 -1.0, NewSafety, substep_no-1, 0);
>> 949 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
>> 950 -1.0, NewSafety, substep_no, 0);
>> 951 G4cout << G4endl;
>> 952
>> 953 #ifdef FUTURE_CORRECTION
>> 954 // Attempt to correct the results of the method // FIX - TODO
>> 955
>> 956 if ( ! found_approximate_intersection )
>> 957 {
>> 958 recalculatedEndPoint = true;
>> 959 // Return the further valid intersection point -- potentially A ??
>> 960 // JA/19 Jan 2006
>> 961 IntersectedOrRecalculatedFT = CurrentA_PointVelocity;
>> 962
>> 963 G4cout << "WARNING - G4PropagatorInField::LocateIntersectionPoint()"
>> 964 << G4endl
>> 965 << " Did not convergence after " << substep_no
>> 966 << " substeps." << G4endl;
>> 967 G4cout << " The endpoint was adjused to pointA resulting"
>> 968 << G4endl
>> 969 << " from the last substep: " << CurrentA_PointVelocity
>> 970 << G4endl;
>> 971 }
>> 972 #endif
>> 973
>> 974 G4cout.precision( 10 );
>> 975 G4double done_len = CurrentA_PointVelocity.GetCurveLength();
>> 976 G4double full_len = CurveEndPointVelocity.GetCurveLength();
>> 977 G4cout << "ERROR - G4PropagatorInField::LocateIntersectionPoint()"
>> 978 << G4endl
>> 979 << " Undertaken only length: " << done_len
>> 980 << " out of " << full_len << " required." << G4endl;
>> 981 G4cout << " Remaining length = " << full_len - done_len << G4endl;
>> 982
>> 983 G4Exception("G4PropagatorInField::LocateIntersectionPoint()",
>> 984 "UnableToLocateIntersection", FatalException,
>> 985 "Too many substeps while trying to locate intersection.");
>> 986 }
>> 987 else if( substep_no >= warn_substeps )
>> 988 {
>> 989 int oldprc= G4cout.precision( 10 );
>> 990 G4cout << "WARNING - G4PropagatorInField::LocateIntersectionPoint()"
>> 991 << G4endl
>> 992 << " Undertaken length: "
>> 993 << CurrentB_PointVelocity.GetCurveLength();
>> 994 G4cout << " - Needed: " << substep_no << " substeps." << G4endl
>> 995 << " Warning level = " << warn_substeps
>> 996 << " and maximum substeps = " << max_substeps << G4endl;
>> 997 G4Exception("G4PropagatorInField::LocateIntersectionPoint()",
>> 998 "DifficultyToLocateIntersection", JustWarning,
>> 999 "Many substeps while trying to locate intersection.");
>> 1000 G4cout.precision( oldprc );
>> 1001 }
>> 1002
>> 1003 return !there_is_no_intersection; // Success or failure
>> 1004 }
>> 1005
>> 1006 ///////////////////////////////////////////////////////////////////////////
531 // 1007 //
>> 1008 // Dumps status of propagator.
>> 1009
532 void 1010 void
533 G4PropagatorInField::printStatus( const G4Fiel << 1011 G4PropagatorInField::printStatus( const G4FieldTrack& StartFT,
534 const G4Fiel << 1012 const G4FieldTrack& CurrentFT,
535 G4doub << 1013 G4double requestStep,
536 G4doub << 1014 G4double safety,
537 G4int << 1015 G4int stepNo,
538 G4VPhy << 1016 G4VPhysicalVolume* startVolume)
539 { 1017 {
540 const G4int verboseLevel = fVerboseLevel; << 1018 const G4int verboseLevel= fVerboseLevel;
541 const G4ThreeVector StartPosition = St 1019 const G4ThreeVector StartPosition = StartFT.GetPosition();
542 const G4ThreeVector StartUnitVelocity = St 1020 const G4ThreeVector StartUnitVelocity = StartFT.GetMomentumDir();
543 const G4ThreeVector CurrentPosition = Cu 1021 const G4ThreeVector CurrentPosition = CurrentFT.GetPosition();
544 const G4ThreeVector CurrentUnitVelocity = Cu 1022 const G4ThreeVector CurrentUnitVelocity = CurrentFT.GetMomentumDir();
545 1023
546 G4double step_len = CurrentFT.GetCurveLength 1024 G4double step_len = CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
547 <<
548 G4long oldprec; // cout/cerr precision set <<
549 1025
550 if( ((stepNo == 0) && (verboseLevel <3)) || << 1026 if( ((stepNo == 0) && (verboseLevel <3))
>> 1027 || (verboseLevel >= 3) )
551 { 1028 {
552 oldprec = G4cout.precision(4); << 1029 static G4int noPrecision= 4;
>> 1030 G4cout.precision(noPrecision);
>> 1031 // G4cout.setf(ios_base::fixed,ios_base::floatfield);
>> 1032 G4cout << std::setw( 6) << " "
>> 1033 << std::setw( 25) << " Current Position and Direction" << " "
>> 1034 << G4endl;
553 G4cout << std::setw( 5) << "Step#" 1035 G4cout << std::setw( 5) << "Step#"
554 << std::setw(10) << " s " << " " 1036 << std::setw(10) << " s " << " "
555 << std::setw(10) << "X(mm)" << " " 1037 << std::setw(10) << "X(mm)" << " "
556 << std::setw(10) << "Y(mm)" << " " 1038 << std::setw(10) << "Y(mm)" << " "
557 << std::setw(10) << "Z(mm)" << " " 1039 << std::setw(10) << "Z(mm)" << " "
558 << std::setw( 7) << " N_x " << " " 1040 << std::setw( 7) << " N_x " << " "
559 << std::setw( 7) << " N_y " << " " 1041 << std::setw( 7) << " N_y " << " "
560 << std::setw( 7) << " N_z " << " " 1042 << std::setw( 7) << " N_z " << " " ;
561 G4cout << std::setw( 7) << " Delta|N|" << << 1043 // << G4endl;
>> 1044 G4cout // << " >>> "
>> 1045 << std::setw( 7) << " Delta|N|" << " "
>> 1046 // << std::setw( 7) << " Delta(N_z) " << " "
562 << std::setw( 9) << "StepLen" << " 1047 << std::setw( 9) << "StepLen" << " "
563 << std::setw(12) << "StartSafety" < 1048 << std::setw(12) << "StartSafety" << " "
564 << std::setw( 9) << "PhsStep" << " 1049 << std::setw( 9) << "PhsStep" << " ";
565 if( startVolume != nullptr ) << 1050 if( startVolume ) {
566 { G4cout << std::setw(18) << "NextVolume << 1051 G4cout << std::setw(18) << "NextVolume" << " ";
567 G4cout.precision(oldprec); << 1052 }
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; 1053 G4cout << G4endl;
605 } 1054 }
606 else // if( verboseLevel > 3 ) << 1055 if((stepNo == 0) && (verboseLevel <=3)){
607 { << 1056 // Recurse to print the start values
608 // Multi-line output << 1057 //
609 << 1058 printStatus( StartFT, StartFT, -1.0, safety, -1, startVolume);
610 G4cout << "Step taken was " << step_len << 1059 }
611 << " out of PhysicalStep = " << re << 1060 if( verboseLevel <= 3 )
612 G4cout << "Final safety is: " << safety << << 1061 {
613 G4cout << "Chord length = " << (CurrentPos << 1062 if( stepNo >= 0)
614 << G4endl; << 1063 G4cout << std::setw( 4) << stepNo << " ";
615 G4cout << G4endl; << 1064 else
616 } << 1065 G4cout << std::setw( 5) << "Start" ;
>> 1066 G4cout.precision(8);
>> 1067 G4cout << std::setw(10) << CurrentFT.GetCurveLength() << " ";
>> 1068 G4cout.precision(8);
>> 1069 G4cout << std::setw(10) << CurrentPosition.x() << " "
>> 1070 << std::setw(10) << CurrentPosition.y() << " "
>> 1071 << std::setw(10) << CurrentPosition.z() << " ";
>> 1072 G4cout.precision(4);
>> 1073 G4cout << std::setw( 7) << CurrentUnitVelocity.x() << " "
>> 1074 << std::setw( 7) << CurrentUnitVelocity.y() << " "
>> 1075 << std::setw( 7) << CurrentUnitVelocity.z() << " ";
>> 1076 // G4cout << G4endl;
>> 1077 // G4cout << " >>> " ;
>> 1078 G4cout.precision(3);
>> 1079 G4cout << std::setw( 7) << CurrentFT.GetMomentum().mag()- StartFT.GetMomentum().mag() << " ";
>> 1080 // << std::setw( 7) << CurrentUnitVelocity.z() - InitialUnitVelocity.z() << " ";
>> 1081 G4cout << std::setw( 9) << step_len << " ";
>> 1082 G4cout << std::setw(12) << safety << " ";
>> 1083 if( requestStep != -1.0 )
>> 1084 G4cout << std::setw( 9) << requestStep << " ";
>> 1085 else
>> 1086 G4cout << std::setw( 9) << "Init/NotKnown" << " ";
>> 1087
>> 1088 if( startVolume != 0)
>> 1089 {
>> 1090 G4cout << std::setw(12) << startVolume->GetName() << " ";
>> 1091 }
>> 1092 #if 0
>> 1093 else
>> 1094 {
>> 1095 if( step_len != -1 )
>> 1096 G4cout << std::setw(12) << "OutOfWorld" << " ";
>> 1097 else
>> 1098 G4cout << std::setw(12) << "NotGiven" << " ";
>> 1099 }
>> 1100 #endif
>> 1101
>> 1102 G4cout << G4endl;
>> 1103 }
>> 1104 else // if( verboseLevel > 3 )
>> 1105 {
>> 1106 // Multi-line output
>> 1107
>> 1108 G4cout << "Step taken was " << step_len
>> 1109 << " out of PhysicalStep= " << requestStep << G4endl;
>> 1110 G4cout << "Final safety is: " << safety << G4endl;
>> 1111
>> 1112 G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
>> 1113 << G4endl;
>> 1114 G4cout << G4endl;
>> 1115 }
617 } 1116 }
618 1117
619 // ------------------------------------------- << 1118 ///////////////////////////////////////////////////////////////////////////
620 // Prints Step diagnostics <<
621 // 1119 //
>> 1120 // Prints Step diagnostics
>> 1121
622 void 1122 void
623 G4PropagatorInField::PrintStepLengthDiagnostic 1123 G4PropagatorInField::PrintStepLengthDiagnostic(
624 G4double CurrentProp 1124 G4double CurrentProposedStepLength,
625 G4double decreaseFac 1125 G4double decreaseFactor,
626 G4double stepTrial, 1126 G4double stepTrial,
627 const G4FieldTrack& ) 1127 const G4FieldTrack& )
628 { 1128 {
629 G4long iprec= G4cout.precision(8); << 1129 G4cout << " PiF: NoZeroStep= " << fNoZeroStep
630 G4cout << " " << std::setw(12) << " PiF: NoZ << 1130 << " CurrentProposedStepLength= " << CurrentProposedStepLength
631 << " " << std::setw(20) << " CurrentP << 1131 << " Full_curvelen_last=" << fFull_CurveLen_of_LastAttempt
632 << " " << std::setw(18) << " Full_cur << 1132 << " last proposed step-length= " << fLast_ProposedStepLength
633 << " " << std::setw(18) << " last pro << 1133 << " decreate factor = " << decreaseFactor
634 << " " << std::setw(18) << " decrease << 1134 << " step trial = " << stepTrial
635 << " " << std::setw(15) << " step tri <<
636 << G4endl; 1135 << G4endl;
>> 1136 }
637 1137
638 G4cout << " " << std::setw(10) << fNoZeroSte << 1138 G4bool
639 << " " << std::setw(20) << CurrentPro << 1139 G4PropagatorInField::IntersectChord( G4ThreeVector StartPointA,
640 << " " << std::setw(18) << fFull_Curv << 1140 G4ThreeVector EndPointB,
641 << " " << std::setw(18) << fLast_Prop << 1141 G4double &NewSafety,
642 << " " << std::setw(18) << decreaseFa << 1142 G4double &LinearStepLength,
643 << " " << std::setw(15) << stepTrial << 1143 G4ThreeVector &IntersectionPoint
644 << G4endl; << 1144 )
645 G4cout.precision( iprec ); << 1145 {
>> 1146 // Calculate the direction and length of the chord AB
>> 1147 G4ThreeVector ChordAB_Vector = EndPointB - StartPointA;
>> 1148 G4double ChordAB_Length = ChordAB_Vector.mag(); // Magnitude (norm)
>> 1149 G4ThreeVector ChordAB_Dir = ChordAB_Vector.unit();
>> 1150 G4bool intersects;
>> 1151
>> 1152 G4ThreeVector OriginShift = StartPointA - fPreviousSftOrigin ;
>> 1153 G4double MagSqShift = OriginShift.mag2() ;
>> 1154 G4double currentSafety;
>> 1155 G4bool doCallNav= false;
>> 1156
>> 1157 if( MagSqShift >= sqr(fPreviousSafety) )
>> 1158 {
>> 1159 currentSafety = 0.0 ;
>> 1160 }else{
>> 1161 currentSafety = fPreviousSafety - std::sqrt(MagSqShift) ;
>> 1162 }
>> 1163
>> 1164 if( fUseSafetyForOptimisation && (ChordAB_Length <= currentSafety) )
>> 1165 {
>> 1166 // The Step is guaranteed to be taken
>> 1167
>> 1168 LinearStepLength = ChordAB_Length;
>> 1169 intersects = false;
>> 1170
>> 1171 NewSafety= currentSafety;
>> 1172
>> 1173 #if 0
>> 1174 G4cout << " G4PropagatorInField does not call Navigator::ComputeStep " << G4endl ;
>> 1175 G4cout << " step= " << LinearStepLength << " safety= " << NewSafety << G4endl;
>> 1176 G4cout << " safety: Origin = " << fPreviousSftOrigin << " val= " << fPreviousSafety << G4endl;
>> 1177 #endif
>> 1178 }
>> 1179 else
>> 1180 {
>> 1181 doCallNav= true;
>> 1182 // Check whether any volumes are encountered by the chord AB
>> 1183
>> 1184 // G4cout << " G4PropagatorInField calling Navigator::ComputeStep " << G4endl ;
>> 1185
>> 1186 LinearStepLength =
>> 1187 fNavigator->ComputeStep( StartPointA, ChordAB_Dir,
>> 1188 ChordAB_Length, NewSafety );
>> 1189 intersects = (LinearStepLength <= ChordAB_Length);
>> 1190 // G4Navigator contracts to return k_infinity if len==asked
>> 1191 // and it did not find a surface boundary at that length
>> 1192 LinearStepLength = std::min( LinearStepLength, ChordAB_Length);
>> 1193
>> 1194 // G4cout << " G4PiF got step= " << LinearStepLength << " safety= " << NewSafety << G4endl;
>> 1195
>> 1196 // Save the last calculated safety!
>> 1197 fPreviousSftOrigin = StartPointA;
>> 1198 fPreviousSafety= NewSafety;
>> 1199
>> 1200 if( intersects ){
>> 1201 // Intersection Point of chord AB and either volume A's surface
>> 1202 // or a daughter volume's surface ..
>> 1203 IntersectionPoint = StartPointA + LinearStepLength * ChordAB_Dir;
>> 1204 }
>> 1205 }
>> 1206
>> 1207 #ifdef DEBUG_INTERSECTS_CHORD
>> 1208 // printIntersection(
>> 1209 // StartPointA, EndPointB, LinearStepLength, IntersectionPoint, NewSafety
>> 1210
>> 1211 G4cout << " G4PropagatorInField::IntersectChord reports " << G4endl;
>> 1212 G4cout << " PiF-IC> "
>> 1213 << "Start=" << std::setw(12) << StartPointA << " "
>> 1214 << "End= " << std::setw(8) << EndPointB << " "
>> 1215 << "StepIn=" << std::setw(8) << LinearStepLength << " "
>> 1216 << "NewSft=" << std::setw(8) << NewSafety << " "
>> 1217 << "CallNav=" << doCallNav << " "
>> 1218 << "Intersects " << intersects << " ";
>> 1219 if( intersects )
>> 1220 G4cout << "IntrPt=" << std::setw(8) << IntersectionPoint << " " ;
>> 1221 G4cout << G4endl;
>> 1222 #endif
>> 1223
>> 1224 return intersects;
>> 1225 }
>> 1226
>> 1227 // --------------------- oooo000000000000oooo ----------------------------
>> 1228
>> 1229 G4FieldTrack G4PropagatorInField::
>> 1230 ReEstimateEndpoint( const G4FieldTrack &CurrentStateA,
>> 1231 const G4FieldTrack &EstimatedEndStateB,
>> 1232 G4double linearDistSq,
>> 1233 G4double curveDist
>> 1234 )
>> 1235 {
>> 1236 // G4double checkCurveDist= EstimatedEndStateB.GetCurveLength()
>> 1237 // - CurrentStateA.GetCurveLength();
>> 1238 // G4double checkLinDistSq= (EstimatedEndStateB.GetPosition()
>> 1239 // - CurrentStateA.GetPosition() ).mag2();
>> 1240
>> 1241 G4FieldTrack newEndPoint( CurrentStateA );
>> 1242 G4MagInt_Driver* integrDriver= GetChordFinder()->GetIntegrationDriver();
>> 1243
>> 1244 G4FieldTrack retEndPoint( CurrentStateA );
>> 1245 G4bool goodAdvance;
>> 1246 G4int itrial=0;
>> 1247 const G4int no_trials= 20;
>> 1248
>> 1249 G4double endCurveLen= EstimatedEndStateB.GetCurveLength();
>> 1250 do
>> 1251 {
>> 1252 G4double currentCurveLen= newEndPoint.GetCurveLength();
>> 1253 G4double advanceLength= endCurveLen - currentCurveLen ;
>> 1254 if (std::abs(advanceLength)<kCarTolerance)
>> 1255 {
>> 1256 advanceLength=(EstimatedEndStateB.GetPosition()
>> 1257 -newEndPoint.GetPosition()).mag();
>> 1258 }
>> 1259 goodAdvance=
>> 1260 integrDriver->AccurateAdvance(newEndPoint, advanceLength, fEpsilonStep);
>> 1261 // ***************
>> 1262 }
>> 1263 while( !goodAdvance && (++itrial < no_trials) );
>> 1264
>> 1265 if( goodAdvance )
>> 1266 {
>> 1267 retEndPoint= newEndPoint;
>> 1268 }
>> 1269 else
>> 1270 {
>> 1271 retEndPoint= EstimatedEndStateB; // Could not improve without major work !!
>> 1272 }
>> 1273
>> 1274 // All the work is done
>> 1275 // below are some diagnostics only -- before the return!
>> 1276 //
>> 1277 static const G4String MethodName("G4PropagatorInField::ReEstimateEndpoint");
>> 1278
>> 1279 #ifdef G4VERBOSE
>> 1280 G4int latest_good_trials=0;
>> 1281 if( itrial > 1)
>> 1282 {
>> 1283 if( fVerboseLevel > 0 )
>> 1284 {
>> 1285 G4cout << MethodName << " called - goodAdv= " << goodAdvance
>> 1286 << " trials = " << itrial
>> 1287 << " previous good= " << latest_good_trials
>> 1288 << G4endl;
>> 1289 }
>> 1290 latest_good_trials=0;
>> 1291 }
>> 1292 else
>> 1293 {
>> 1294 latest_good_trials++;
>> 1295 }
>> 1296 #endif
>> 1297
>> 1298 #ifdef G4DEBUG_FIELD
>> 1299 G4double lengthDone = newEndPoint.GetCurveLength()
>> 1300 - CurrentStateA.GetCurveLength();
>> 1301 if( !goodAdvance )
>> 1302 {
>> 1303 if( fVerboseLevel >= 3 )
>> 1304 {
>> 1305 G4cout << MethodName << "> AccurateAdvance failed " ;
>> 1306 G4cout << " in " << itrial << " integration trials/steps. " << G4endl;
>> 1307 G4cout << " It went only " << lengthDone << " instead of " << curveDist
>> 1308 << " -- a difference of " << curveDist - lengthDone << G4endl;
>> 1309 G4cout << " ReEstimateEndpoint> Reset endPoint to original value!"
>> 1310 << G4endl;
>> 1311 }
>> 1312 }
>> 1313
>> 1314 static G4int noInaccuracyWarnings = 0;
>> 1315 G4int maxNoWarnings = 10;
>> 1316 if ( (noInaccuracyWarnings < maxNoWarnings )
>> 1317 || (fVerboseLevel > 1) )
>> 1318 {
>> 1319 G4cerr << "G4PropagatorInField::LocateIntersectionPoint():"
>> 1320 << G4endl
>> 1321 << " Warning: Integration inaccuracy requires"
>> 1322 << " an adjustment in the step's endpoint." << G4endl
>> 1323 << " Two mid-points are further apart than their"
>> 1324 << " curve length difference" << G4endl
>> 1325 << " Dist = " << std::sqrt(linearDistSq)
>> 1326 << " curve length = " << curveDist << G4endl;
>> 1327 G4cerr << " Correction applied is "
>> 1328 << (newEndPoint.GetPosition()-EstimatedEndStateB.GetPosition()).mag()
>> 1329 << G4endl;
>> 1330 }
>> 1331 #else
>> 1332 // Statistics on the RMS value of the corrections
>> 1333
>> 1334 static G4int noCorrections=0;
>> 1335 static G4double sumCorrectionsSq = 0;
>> 1336 noCorrections++;
>> 1337 if( goodAdvance )
>> 1338 {
>> 1339 sumCorrectionsSq += (EstimatedEndStateB.GetPosition() -
>> 1340 newEndPoint.GetPosition()).mag2();
>> 1341 }
>> 1342 linearDistSq -= curveDist; // To use linearDistSq ... !
>> 1343 #endif
>> 1344
>> 1345 return retEndPoint;
646 } 1346 }
647 1347
648 // Access the points which have passed through 1348 // Access the points which have passed through the filter. The
649 // points are stored as ThreeVectors for the i 1349 // points are stored as ThreeVectors for the initial impelmentation
650 // only (jacek 30/10/2002) 1350 // only (jacek 30/10/2002)
651 // Responsibility for deleting the points lies 1351 // Responsibility for deleting the points lies with
652 // SmoothTrajectoryPoint, which is the points' 1352 // SmoothTrajectoryPoint, which is the points' final
653 // destination. The points pointer is set to N 1353 // destination. The points pointer is set to NULL, to ensure that
654 // the points are not re-used in subsequent st 1354 // the points are not re-used in subsequent steps, therefore THIS
655 // METHOD MUST BE CALLED EXACTLY ONCE PER STEP 1355 // METHOD MUST BE CALLED EXACTLY ONCE PER STEP. (jacek 08/11/2002)
656 1356
657 std::vector<G4ThreeVector>* 1357 std::vector<G4ThreeVector>*
658 G4PropagatorInField::GimmeTrajectoryVectorAndF 1358 G4PropagatorInField::GimmeTrajectoryVectorAndForgetIt() const
659 { 1359 {
660 // NB, GimmeThePointsAndForgetThem really fo 1360 // NB, GimmeThePointsAndForgetThem really forgets them, so it can
661 // only be called (exactly) once for each st 1361 // only be called (exactly) once for each step.
662 1362
663 if (fpTrajectoryFilter != nullptr) << 1363 if (fpTrajectoryFilter)
664 { 1364 {
665 return fpTrajectoryFilter->GimmeThePointsA 1365 return fpTrajectoryFilter->GimmeThePointsAndForgetThem();
666 } 1366 }
667 return nullptr; << 1367 else
>> 1368 {
>> 1369 return 0;
>> 1370 }
668 } 1371 }
669 1372
670 // ------------------------------------------- <<
671 // <<
672 void 1373 void
673 G4PropagatorInField::SetTrajectoryFilter(G4VCu 1374 G4PropagatorInField::SetTrajectoryFilter(G4VCurvedTrajectoryFilter* filter)
674 { 1375 {
675 fpTrajectoryFilter = filter; 1376 fpTrajectoryFilter = filter;
676 } 1377 }
677 1378
678 // ------------------------------------------- <<
679 // <<
680 void G4PropagatorInField::ClearPropagatorState 1379 void G4PropagatorInField::ClearPropagatorState()
681 { 1380 {
682 // Goal: Clear all memory of previous steps, 1381 // Goal: Clear all memory of previous steps, cached information
683 1382
684 fParticleIsLooping = false; << 1383 fParticleIsLooping= false;
685 fNoZeroStep = 0; << 1384 fNoZeroStep= 0;
686 1385
687 fSetFieldMgr = false; // Has field-manager <<
688 fEpsilonStep= 1.0e-5; // Relative accuracy <<
689 <<
690 End_PointAndTangent= G4FieldTrack( G4ThreeVe 1386 End_PointAndTangent= G4FieldTrack( G4ThreeVector(0.,0.,0.),
691 G4ThreeVe 1387 G4ThreeVector(0.,0.,0.),
692 0.0,0.0,0 1388 0.0,0.0,0.0,0.0,0.0);
693 fFull_CurveLen_of_LastAttempt = -1; 1389 fFull_CurveLen_of_LastAttempt = -1;
694 fLast_ProposedStepLength = -1; 1390 fLast_ProposedStepLength = -1;
695 1391
696 fPreviousSftOrigin= G4ThreeVector(0.,0.,0.); 1392 fPreviousSftOrigin= G4ThreeVector(0.,0.,0.);
697 fPreviousSafety= 0.0; 1393 fPreviousSafety= 0.0;
698 <<
699 fNewTrack = true; <<
700 } 1394 }
701 1395
702 // ------------------------------------------- <<
703 // <<
704 G4FieldManager* G4PropagatorInField:: 1396 G4FieldManager* G4PropagatorInField::
705 FindAndSetFieldManager( G4VPhysicalVolume* pCu << 1397 FindAndSetFieldManager( G4VPhysicalVolume* pCurrentPhysicalVolume)
706 { 1398 {
707 G4FieldManager* currentFieldMgr; 1399 G4FieldManager* currentFieldMgr;
708 1400
709 currentFieldMgr = fDetectorFieldMgr; 1401 currentFieldMgr = fDetectorFieldMgr;
710 if( pCurrentPhysicalVolume != nullptr ) << 1402 if( pCurrentPhysicalVolume)
711 { 1403 {
712 G4FieldManager *pRegionFieldMgr = nullptr << 1404 G4FieldManager *newFieldMgr = 0;
713 G4LogicalVolume* pLogicalVol = pCurrentPh << 1405 newFieldMgr= pCurrentPhysicalVolume->GetLogicalVolume()->GetFieldManager();
714 << 1406 if ( newFieldMgr )
715 if( pLogicalVol != nullptr ) << 1407 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 } 1408 }
738 fCurrentFieldMgr = currentFieldMgr; << 1409 fCurrentFieldMgr= currentFieldMgr;
739 1410
740 // Flag that field manager has been set << 1411 // Flag that field manager has been set.
741 // << 1412 fSetFieldMgr= true;
742 fSetFieldMgr = true; <<
743 1413
744 return currentFieldMgr; 1414 return currentFieldMgr;
745 } 1415 }
746 1416
747 // ------------------------------------------- <<
748 // <<
749 G4int G4PropagatorInField::SetVerboseLevel( G4 1417 G4int G4PropagatorInField::SetVerboseLevel( G4int level )
750 { 1418 {
751 G4int oldval = fVerboseLevel; << 1419 G4int oldval= fVerboseLevel;
752 fVerboseLevel = level; << 1420 fVerboseLevel= level;
753 1421
754 // Forward the verbose level 'reduced' to Ch 1422 // Forward the verbose level 'reduced' to ChordFinder,
755 // MagIntegratorDriver ... ? 1423 // MagIntegratorDriver ... ?
756 // 1424 //
757 auto integrDriver = GetChordFinder()->GetInt << 1425 G4MagInt_Driver* integrDriver= GetChordFinder()->GetIntegrationDriver();
758 integrDriver->SetVerboseLevel( fVerboseLevel << 1426 integrDriver->SetVerboseLevel( fVerboseLevel - 2 );
759 G4cout << "Set Driver verbosity to " << fVer 1427 G4cout << "Set Driver verbosity to " << fVerboseLevel - 2 << G4endl;
760 1428
761 return oldval; 1429 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 } 1430 }
886 1431