| 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 <<
27 // 26 //
28 // class description: << 27 // $Id: G4PropagatorInField.hh,v 1.9 2006/06/29 18:36:15 gunter Exp $
>> 28 // GEANT4 tag $Name: geant4-08-01-patch-01 $
>> 29 //
>> 30 // class G4PropagatorInField
>> 31 //
>> 32 // Class description:
29 // 33 //
30 // This class performs the navigation/propagat 34 // This class performs the navigation/propagation of a particle/track
31 // in a magnetic field. The field is in genera 35 // in a magnetic field. The field is in general non-uniform.
32 // For the calculation of the path, it relies 36 // For the calculation of the path, it relies on the class G4ChordFinder.
33 << 37 //
>> 38 // Key Method:
>> 39 // ComputeStep(..)
34 // History: 40 // History:
35 // ------- 41 // -------
36 // 25.10.96 John Apostolakis, design and impl 42 // 25.10.96 John Apostolakis, design and implementation
37 // 25.03.97 John Apostolakis, adaptation for 43 // 25.03.97 John Apostolakis, adaptation for G4Transportation and cleanup
38 // 8.11.02 John Apostolakis, changes to enab 44 // 8.11.02 John Apostolakis, changes to enable use of safety in intersecting
39 // ------------------------------------------- 45 // ---------------------------------------------------------------------------
>> 46
40 #ifndef G4PropagatorInField_hh 47 #ifndef G4PropagatorInField_hh
41 #define G4PropagatorInField_hh 1 48 #define G4PropagatorInField_hh 1
42 49
43 #include "G4Types.hh" 50 #include "G4Types.hh"
44 51
45 #include <vector> 52 #include <vector>
46 53
47 #include "G4FieldTrack.hh" 54 #include "G4FieldTrack.hh"
48 #include "G4FieldManager.hh" 55 #include "G4FieldManager.hh"
49 #include "G4VIntersectionLocator.hh" <<
50 <<
51 class G4ChordFinder; 56 class G4ChordFinder;
52 57
53 class G4Navigator; 58 class G4Navigator;
54 class G4VPhysicalVolume; 59 class G4VPhysicalVolume;
55 class G4VCurvedTrajectoryFilter; 60 class G4VCurvedTrajectoryFilter;
56 61
57 class G4PropagatorInField 62 class G4PropagatorInField
58 { 63 {
59 64
60 public: // with description 65 public: // with description
61 66
62 G4PropagatorInField( G4Navigator* theNaviga << 67 G4PropagatorInField( G4Navigator *theNavigator,
63 G4FieldManager* detect << 68 G4FieldManager *detectorFieldMgr );
64 G4VIntersectionLocator <<
65 ~G4PropagatorInField(); 69 ~G4PropagatorInField();
66 70
67 G4double ComputeStep( G4FieldTrack& pFieldT << 71 G4double ComputeStep( G4FieldTrack &pFieldTrack,
68 G4double pCurrentProp << 72 G4double pCurrentProposedStepLength,
69 G4double& pNewSafety, << 73 G4double &pNewSafety,
70 G4VPhysicalVolume* pP << 74 G4VPhysicalVolume *pPhysVol=0 );
71 G4bool canRelaxDeltaC <<
72 // Compute the next geometric Step 75 // Compute the next geometric Step
73 76
74 inline G4ThreeVector EndPosition() const; << 77 inline G4ThreeVector EndPosition() const;
75 inline G4ThreeVector EndMomentumDir() const << 78 inline G4ThreeVector EndMomentumDir() const;
76 inline G4bool IsParticleLooping() co << 79 inline G4bool IsParticleLooping() const;
77 // Return the state after the Step 80 // Return the state after the Step
78 81
79 inline G4double GetEpsilonStep() const; << 82 inline G4double GetEpsilonStep() const;
80 // Relative accuracy for current Step (Ca 83 // Relative accuracy for current Step (Calc.)
81 inline void SetEpsilonStep(G4double new << 84 inline void SetEpsilonStep(G4double newEps);
82 // The ratio DeltaOneStep()/h_current_ste 85 // The ratio DeltaOneStep()/h_current_step
83 86
84 G4FieldManager* FindAndSetFieldManager(G4VP << 87 inline void SetChargeMomentumMass( G4double charge, // in e+ units
>> 88 G4double momentum, // in Geant4 units
>> 89 G4double pMass );
>> 90 // Inform this and all associated classes of q, p, m
>> 91
>> 92 G4FieldManager* FindAndSetFieldManager(G4VPhysicalVolume* pCurrentPhysVol);
85 // Set (and return) the correct field man 93 // Set (and return) the correct field manager (global or local),
86 // if it exists. << 94 // if it exists.
87 // Should be called before ComputeStep is 95 // Should be called before ComputeStep is called;
88 // Currently, ComputeStep will call it, i << 96 // - currently, ComputeStep will call it, if it has not been called.
89 97
90 inline G4ChordFinder* GetChordFinder(); 98 inline G4ChordFinder* GetChordFinder();
91 99
92 G4int SetVerboseLevel( G4int verbose << 100 G4int SetVerboseLevel( G4int verbose );
93 inline G4int GetVerboseLevel() const; << 101 inline G4int GetVerboseLevel() const;
94 inline G4int Verbose() const; << 102 inline G4int Verbose() const;
95 inline void CheckMode(G4bool mode); << 103
96 << 104 inline G4int GetMaxLoopCount() const;
97 inline void SetVerboseTrace( G4bool enabl << 105 inline void SetMaxLoopCount( G4int new_max );
98 inline G4bool GetVerboseTrace(); << 106 // A maximum for the number of steps that a (looping) particle can take.
99 // Tracing key parts of Compute Step << 107
100 << 108 void printStatus( const G4FieldTrack& startFT,
101 inline G4int GetMaxLoopCount() const; << 109 const G4FieldTrack& currentFT,
102 inline void SetMaxLoopCount( G4int new_max << 110 G4double requestStep,
103 // A maximum for the number of substeps t << 111 G4double safety,
104 // Above this number it is signaled as << 112 G4int step,
105 << 113 G4VPhysicalVolume* startVolume);
106 void printStatus( const G4FieldTrack& <<
107 const G4FieldTrack& <<
108 G4double <<
109 G4double <<
110 G4int <<
111 G4VPhysicalVolume* <<
112 // Print Method - useful mostly for debug 114 // Print Method - useful mostly for debugging.
113 115
114 inline G4FieldTrack GetEndState() const; 116 inline G4FieldTrack GetEndState() const;
115 117
116 inline G4double GetMinimumEpsilonStep() con << 118 // The following methods are now obsolescent but *for now* will work
117 inline void SetMinimumEpsilonStep( G4do << 119 // They are being replaced by same-name methods in G4FieldManager,
118 inline G4double GetMaximumEpsilonStep() con << 120 // allowing the specialisation in different volumes.
119 inline void SetMaximumEpsilonStep( G4do << 121 // Their new behaviour is to change the values for the global field manager.
120 // The 4 above methods are now obsolescen << 122 inline G4double GetMinimumEpsilonStep() const;
121 // They are being replaced by same-name m << 123 inline void SetMinimumEpsilonStep( G4double newEpsMin );
122 // allowing the specialisation in differe << 124 // Minimum for Relative accuracy of any Step
123 // Their new behaviour is to change the v <<
124 // manager <<
125 <<
126 void SetLargestAcceptableStep( G4double <<
127 G4double GetLargestAcceptableStep(); <<
128 void ResetLargestAcceptableStep(); <<
129 // Obtain / change the size of the larges <<
130 // Reset method uses the world volume's <<
131 <<
132 G4double GetMaxStepSizeMultiplier(); <<
133 void SetMaxStepSizeMultiplier(G4double <<
134 // Control extra Multiplier parameter for <<
135 G4double GetMinBigDistance(); <<
136 void SetMinBigDistance(G4double val); <<
137 // Control minimum 'directional' distance <<
138 <<
139 void SetTrajectoryFilter(G4VCurvedTrajector <<
140 // Set the filter that examines & stores <<
141 // curved trajectory points. Currently o <<
142 <<
143 std::vector<G4ThreeVector>* GimmeTrajectory <<
144 // Access the points which have passed by <<
145 // Responsibility for deleting the points <<
146 // This method MUST BE called exactly ONC <<
147 <<
148 void ClearPropagatorState(); <<
149 // Clear all the State of this class and <<
150 // --> the current field manager & chord <<
151 125
152 inline void SetDetectorFieldManager( G4Fiel << 126 inline G4double GetMaximumEpsilonStep() const;
153 // Update this (dangerous) state -- for t << 127 inline void SetMaximumEpsilonStep( G4double newEpsMax );
154 <<
155 inline void SetUseSafetyForOptimization( <<
156 inline G4bool GetUseSafetyForOptimization() <<
157 // Toggle & view parameter for using safe <<
158 // unneccesary calls to navigator (thus ' <<
159 inline G4bool IntersectChord( const G4Three <<
160 const G4Three <<
161 G4doubl <<
162 G4doubl <<
163 G4Three <<
164 // Intersect the chord from StartPointA t <<
165 // and return whether an intersection occ <<
166 // NOTE: Safety is changed! <<
167 128
168 inline G4bool IsFirstStepInVolume(); << 129 inline void SetLargestAcceptableStep( G4double newBigDist );
169 inline G4bool IsLastStepInVolume(); << 130 inline G4double GetLargestAcceptableStep();
170 inline void PrepareNewTrack(); <<
171 <<
172 inline G4VIntersectionLocator* GetIntersect <<
173 inline void SetIntersectionLocator(G4VInter <<
174 // Change or get the object which calcula <<
175 // intersection point with the next bound <<
176 <<
177 inline G4int GetIterationsToIncreaseChordDi <<
178 inline void SetIterationsToIncreaseChordDi <<
179 // Control the parameter which enables th <<
180 // which ensures that chord segments ar <<
181 // their sagitta is small than delta-ch <<
182 // The Set method increases the value of <<
183 // doubling it once the number of itera <<
184 // value of 'IncreaseChordDistanceThres <<
185 // again every time the iteration count <<
186 // value. <<
187 // Note: delta-chord is reset to its orig <<
188 // each call to ComputeStep. <<
189 131
190 public: // without description << 132 public: // with description
>> 133
>> 134 // The following methods are obsolete and will not work --
>> 135 // as they have been replaced by the same methods in G4FieldManager
>> 136 // since Geant4 4.0
>> 137 inline G4double GetDeltaIntersection() const;
>> 138 inline G4double GetDeltaOneStep() const;
>> 139 inline void SetAccuraciesWithDeltaOneStep( G4double deltaOneStep );
>> 140 inline void SetDeltaIntersection( G4double deltaIntersection );
>> 141 inline void SetDeltaOneStep( G4double deltaOneStep );
191 142
192 inline G4double GetDeltaIntersection() cons << 143 public: // without description
193 inline G4double GetDeltaOneStep() const; <<
194 144
195 inline G4FieldManager* GetCurrentFieldManag << 145 inline G4FieldManager* GetCurrentFieldManager();
196 inline G4EquationOfMotion* GetCurrentEquati <<
197 // Auxiliary methods - their results can <<
198 146
199 inline void SetNavigatorForPropagating(G4Na << 147 public: // no description
200 inline G4Navigator* GetNavigatorForPropagat <<
201 148
202 inline void SetThresholdNoZeroStep( G4int n 149 inline void SetThresholdNoZeroStep( G4int noAct,
203 G4int n 150 G4int noHarsh,
204 G4int n 151 G4int noAbandon );
205 inline G4int GetThresholdNoZeroSteps( G4int 152 inline G4int GetThresholdNoZeroSteps( G4int i );
206 153
207 inline G4double GetZeroStepThreshold(); << 154 public: // with description
208 inline void SetZeroStepThreshold( G4dou << 155 //
209 << 156 void SetTrajectoryFilter(G4VCurvedTrajectoryFilter* filter);
210 void RefreshIntersectionLocator(); << 157 // Set the filter that examines & stores 'intermediate'
211 // Update the Locator with parameters fro << 158 // curved trajectory points. Currently only position is stored.
212 // and from current field manager << 159
>> 160 std::vector<G4ThreeVector>* GimmeTrajectoryVectorAndForgetIt() const;
>> 161 // Access the points which have passed by the filter.
>> 162 // Responsibility for deleting the points lies with the client.
>> 163 // This method MUST BE called exactly ONCE per step.
>> 164
>> 165 void ClearPropagatorState();
>> 166 // Clear all the State of this class and its current associates
>> 167 // --> the current field manager & chord finder will also be called
>> 168
>> 169 inline void SetDetectorFieldManager( G4FieldManager* newGlobalFieldManager );
>> 170 // Update this (dangerous) state -- for the time being
>> 171
>> 172 inline void SetUseSafetyForOptimization( G4bool );
>> 173 inline G4bool GetUseSafetyForOptimization();
>> 174 // Toggle & view parameter for using safety to discard
>> 175 // unneccesary calls to navigator (thus 'optimising' performance)
>> 176
>> 177 protected: // with description
>> 178
>> 179 G4bool LocateIntersectionPoint(
>> 180 const G4FieldTrack& curveStartPointTangent, // A
>> 181 const G4FieldTrack& curveEndPointTangent, // B
>> 182 const G4ThreeVector& trialPoint, // E
>> 183 G4FieldTrack& intersectPointTangent, // Output
>> 184 G4bool& recalculatedEndPoint); // Out:
>> 185
>> 186 // If such an intersection exists, this function
>> 187 // calculate the intersection point of the true path of the particle
>> 188 // with the surface of the current volume (or of one of its daughters).
>> 189 // (Should use lateral displacement as measure of convergence).
>> 190
>> 191 G4bool IntersectChord( G4ThreeVector StartPointA,
>> 192 G4ThreeVector EndPointB,
>> 193 G4double &NewSafety,
>> 194 G4double &LinearStepLength,
>> 195 G4ThreeVector &IntersectionPoint);
>> 196 // Intersect the chord from StartPointA to EndPointB
>> 197 // and return whether an intersection occurred
213 198
214 protected: // without description << 199 G4FieldTrack ReEstimateEndpoint( const G4FieldTrack &CurrentStateA,
>> 200 const G4FieldTrack &EstimtdEndStateB,
>> 201 G4double linearDistSq,
>> 202 G4double curveDist);
>> 203 // Return new estimate for state after curveDist
>> 204 // starting from CurrentStateA, to replace EstimtdEndStateB,
>> 205 // (and report displacement -- if field is compiled verbose.)
215 206
216 void PrintStepLengthDiagnostic( G4double 207 void PrintStepLengthDiagnostic( G4double currentProposedStepLength,
217 G4double 208 G4double decreaseFactor,
218 G4double 209 G4double stepTrial,
219 const G4FieldTrac 210 const G4FieldTrack& aFieldTrack);
220 <<
221 void ReportLoopingParticle( G4int count, G <<
222 G4double stepRe <<
223 const G4ThreeVe <<
224 G4VPhysicalVolu <<
225 void ReportStuckParticle(G4int noZeroSteps, <<
226 G4double lastTried <<
227 <<
228 private: 211 private:
229 212
230 // ---------------------------------------- << 213 // ----------------------------------------------------------------------
231 // DATA Members << 214 // DATA Members
232 // ---------------------------------------- << 215 // ----------------------------------------------------------------------
233 <<
234 // ======================================= <<
235 // INVARIANTS - Must not change during tra <<
236 <<
237 // ** PARAMETERS ----------- <<
238 G4int fMax_loop_count = 1000; <<
239 // Limit for the number of sub-steps take <<
240 G4int fIncreaseChordDistanceThreshold = 100 <<
241 G4bool fUseSafetyForOptimisation = true; <<
242 // (false) is less sensitive to incorrect <<
243 216
244 // Thresholds for identifying "abnormal" c << 217 G4FieldManager *fDetectorFieldMgr;
245 // << 218 // The Field Manager of the whole Detector. (default)
246 G4int fActionThreshold_NoZeroSteps = 2; <<
247 G4int fSevereActionThreshold_NoZeroSteps = <<
248 G4int fAbandonThreshold_NoZeroSteps = 50; <<
249 G4double fZeroStepThreshold = 0.0; <<
250 // Threshold *length* for counting of tin <<
251 <<
252 // Parameters related to handling of very l <<
253 // occur typically in large volumes with <<
254 G4double fLargestAcceptableStep; <<
255 // Maximum size of a step - for optimizat <<
256 G4double fMaxStepSizeMultiplier = 3; <<
257 // Multiplier for directional exit distan <<
258 G4double fMinBigDistance= 100. ; // * CLHEP <<
259 // Minimum distance added to directional <<
260 // ** End of PARAMETERS ----- <<
261 <<
262 G4double kCarTolerance; <<
263 // Geometrical tolerance defining surfa <<
264 <<
265 G4bool fAllocatedLocator; <<
266 <<
267 // --------------------------------------- <<
268 // ** Dependent Objects - to which work is <<
269 <<
270 G4FieldManager* fDetectorFieldMgr; <<
271 // The Field Manager of the whole Dete <<
272 <<
273 G4VIntersectionLocator* fIntersectionLocato <<
274 // Refines candidate intersection <<
275 <<
276 G4VCurvedTrajectoryFilter* fpTrajectoryFilt <<
277 // The filter encapsulates the algorithm <<
278 // intermediate points should be stored i <<
279 // When it is NULL, no intermediate point <<
280 // Else PIF::ComputeStep must submit (all <<
281 // points it calculates, to this filter. <<
282 219
283 G4Navigator* fNavigator; << 220 G4FieldManager *fCurrentFieldMgr;
284 // Set externally - only by tracking / ru << 221 // The Field Manager of the current volume (may be the one above.)
285 // <<
286 // ** End of Dependent Objects ----------- <<
287 222
288 // End of INVARIANTS << 223 G4Navigator *fNavigator;
289 // ======================================= <<
290 224
291 // STATE information 225 // STATE information
292 // ----------------- 226 // -----------------
293 G4FieldManager* fCurrentFieldMgr; << 227
294 // The Field Manager of the current volu << 228 G4double fEpsilonStep;
295 G4bool fSetFieldMgr = false; // Has it bee << 229 // Relative accuracy for current Step (Calc.)
296 << 230
297 // Parameters of current step << 231 G4FieldTrack End_PointAndTangent;
298 G4double fEpsilonStep; // Relati << 232 // End point storage
299 G4FieldTrack End_PointAndTangent; // End po << 233
300 G4bool fParticleIsLooping = false; << 234 G4bool fParticleIsLooping;
301 G4int fNoZeroStep = 0; // Count << 235
302 << 236 G4int fVerboseLevel;
303 // State used for Optimisation << 237 // For debuging purposes
304 G4double fFull_CurveLen_of_LastAttempt = -1 << 238
305 G4double fLast_ProposedStepLength = -1; << 239 // Limit for the number of sub-steps taken in one call to ComputeStep
306 // Previous step information -- for use i << 240 G4int fMax_loop_count;
307 G4ThreeVector fPreviousSftOrigin; << 241
308 G4double fPreviousSafety = 0.0; << 242 // Variables to keep track of "abnormal" case - which causes loop
309 // Last safety origin & value: for optimi << 243 //
310 << 244 G4int fNoZeroStep; // Counter of zeroStep
311 G4int fVerboseLevel = 0; << 245 G4int fActionThreshold_NoZeroSteps; // Threshold: above this - act
312 G4bool fVerbTracePiF = false; << 246 G4int fSevereActionThreshold_NoZeroSteps; // Threshold to act harshly
313 G4bool fCheck = false; << 247 G4int fAbandonThreshold_NoZeroSteps; // Threshold to abandon
314 // For debugging purposes << 248
315 << 249 G4double fFull_CurveLen_of_LastAttempt;
316 G4bool fFirstStepInVolume = true; << 250 G4double fLast_ProposedStepLength;
317 G4bool fLastStepInVolume = true; << 251 G4double fLargestAcceptableStep;
318 G4bool fNewTrack = true; << 252
>> 253 G4double fCharge, fInitialMomentumModulus, fMass;
>> 254
>> 255 // Last safety origin & value: for optimisation
>> 256 G4ThreeVector fPreviousSftOrigin;
>> 257 G4double fPreviousSafety;
>> 258 G4bool fUseSafetyForOptimisation;
>> 259
>> 260 // Flag whether field manager has been set for the current step
>> 261 G4bool fSetFieldMgr;
>> 262
>> 263 private:
>> 264 // The filter encapsulates the algorithm which selects which
>> 265 // intermediate points should be stored in a trajectory.
>> 266 // When it is NULL, no intermediate points will be stored.
>> 267 // Else PIF::ComputeStep must submit (all) intermediate
>> 268 // points it calculates, to this filter. (jacek 04/11/2002)
>> 269 G4VCurvedTrajectoryFilter* fpTrajectoryFilter;
319 }; 270 };
320 271
321 // Inline methods << 272 // ********************************************************************
322 // << 273 // Inline methods.
>> 274 // ********************************************************************
>> 275
323 #include "G4PropagatorInField.icc" 276 #include "G4PropagatorInField.icc"
324 277
325 #endif 278 #endif
>> 279
326 280