ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MagIntegratorDriver.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MagIntegratorDriver.hh
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4MagInt_Driver
27 //
28 // Class description:
29 //
30 // Provides a driver that talks to the Integrator Stepper, and insures that
31 // the error is within acceptable bounds.
32 
33 // V.Grichine, 07.10.1996 - Created
34 // W.Wander, 28.01.1998 - Added ability for low order integrators
35 // J.Apostolakis, 08.11.2001 - Respect minimum step in AccurateAdvance
36 // --------------------------------------------------------------------
37 #ifndef G4MAGINT_DRIVER_HH
38 #define G4MAGINT_DRIVER_HH
39 
40 #include "G4VIntegrationDriver.hh"
41 #include "G4MagIntegratorStepper.hh"
42 #include "G4ChordFinderDelegate.hh"
43 
44 class G4MagInt_Driver : public G4VIntegrationDriver,
45  public G4ChordFinderDelegate<G4MagInt_Driver>
46 {
47  public: // with description
48 
49  G4MagInt_Driver(G4double hminimum,
50  G4MagIntegratorStepper* pItsStepper,
51  G4int numberOfComponents = 6,
52  G4int statisticsVerbosity = 0);
53  virtual ~G4MagInt_Driver() override;
54  // Constructor, destructor.
55 
56  G4MagInt_Driver(const G4MagInt_Driver&) = delete;
57  G4MagInt_Driver& operator=(const G4MagInt_Driver&) = delete;
58 
59  inline virtual G4double AdvanceChordLimited(G4FieldTrack& track,
60  G4double stepMax,
61  G4double epsStep,
62  G4double chordDistance) override;
63 
64  inline virtual void OnStartTracking() override;
65  inline virtual void OnComputeStep() override {};
66  virtual G4bool DoesReIntegrate() override { return true; }
67 
68  virtual G4bool AccurateAdvance(G4FieldTrack& y_current,
69  G4double hstep,
70  G4double eps, // Requested y_err/hstep
71  G4double hinitial = 0.0) override;
72  // Above drivers for integrator (Runge-Kutta) with stepsize control.
73  // Integrates ODE starting values y_current
74  // from current s (s=s0) to s=s0+h with accuracy eps.
75  // On output ystart is replaced by value at end of interval.
76  // The concept is similar to the odeint routine from NRC p.721-722.
77 
78  virtual G4bool QuickAdvance(G4FieldTrack& y_val, // INOUT
79  const G4double dydx[],
80  G4double hstep,
81  G4double& dchord_step,
82  G4double& dyerr) override;
83  // QuickAdvance just tries one Step - it does not ensure accuracy.
84 
85  G4bool QuickAdvance(G4FieldTrack& y_posvel, // INOUT
86  const G4double dydx[],
87  G4double hstep, // IN
88  G4double& dchord_step,
89  G4double& dyerr_pos_sq,
90  G4double& dyerr_mom_rel_sq );
91  // New QuickAdvance that also just tries one Step
92  // (so also does not ensure accuracy)
93  // but does return the errors in position and
94  // momentum (normalised: Delta_Integration(p^2)/(p^2) )
95 
96  inline G4double GetHmin() const;
97  inline G4double Hmin() const; // Obsolete
98  inline G4double GetSafety() const;
99  inline G4double GetPshrnk() const;
100  inline G4double GetPgrow() const;
101  inline G4double GetErrcon() const;
102  virtual void GetDerivatives(const G4FieldTrack& y_curr, // INput
103  G4double dydx[]) const override; // OUTput
104 
105  virtual void GetDerivatives(const G4FieldTrack& track,
106  G4double dydx[],
107  G4double field[]) const override;
108  // Accessors
109 
110  virtual G4EquationOfMotion* GetEquationOfMotion() override;
111  virtual void SetEquationOfMotion(G4EquationOfMotion* equation) override;
112 
113  virtual void RenewStepperAndAdjust(G4MagIntegratorStepper* pItsStepper) override;
114  // Sets a new stepper pItsStepper for this driver. Then it calls
115  // ReSetParameters to reset its parameters accordingly.
116 
117  inline void ReSetParameters(G4double new_safety = 0.9);
118  // i) sets the exponents (pgrow & pshrnk),
119  // using the current Stepper's order,
120  // ii) sets the safety
121  // ii) calculates "errcon" according to the above values.
122 
123  inline void SetSafety(G4double valS);
124  inline void SetPshrnk(G4double valPs);
125  inline void SetPgrow (G4double valPg);
126  inline void SetErrcon(G4double valEc);
127  // When setting safety or pgrow, errcon will be set to a compatible value.
128 
129  inline G4double ComputeAndSetErrcon();
130 
131  virtual const G4MagIntegratorStepper* GetStepper() const override;
132  virtual G4MagIntegratorStepper* GetStepper() override;
133 
134  void OneGoodStep(G4double ystart[], // Like old RKF45step()
135  const G4double dydx[],
136  G4double& x,
137  G4double htry,
138  G4double eps, // memb variables ?
139  G4double& hdid,
140  G4double& hnext ) ;
141  // This takes one Step that is as large as possible while
142  // satisfying the accuracy criterion of:
143  // yerr < eps * |y_end-y_start|
144 
145  virtual G4double ComputeNewStepSize(G4double errMaxNorm, // normalised
146  G4double hstepCurrent) override;
147  // Taking the last step's normalised error, calculate
148  // a step size for the next step.
149  // Do not limit the next step's size within a factor of the
150  // current one.
151 
152  G4double ComputeNewStepSize_WithinLimits(G4double errMaxNorm, // normalised
153  G4double hstepCurrent);
154  // Taking the last step's normalised error, calculate
155  // a step size for the next step.
156  // Limit the next step's size within a range around the current one.
157 
158  inline G4int GetMaxNoSteps() const;
159  inline void SetMaxNoSteps(G4int val);
160  // Modify and Get the Maximum number of Steps that can be
161  // taken for the integration of a single segment -
162  // (i.e. a single call to AccurateAdvance).
163 
164  public: // without description
165 
166  inline void SetHmin(G4double newval);
167  virtual void SetVerboseLevel(G4int newLevel) override;
168  virtual G4int GetVerboseLevel() const override;
169 
170  inline G4double GetSmallestFraction() const;
171  void SetSmallestFraction( G4double val );
172 
173  protected: // without description
174 
175  void WarnSmallStepSize(G4double hnext, G4double hstep,
176  G4double h, G4double xDone,
177  G4int noSteps);
178 
179  void WarnTooManyStep(G4double x1start, G4double x2end, G4double xCurrent);
180  void WarnEndPointTooFar(G4double endPointDist,
181  G4double hStepSize ,
182  G4double epsilonRelative,
183  G4int debugFlag);
184  // Issue warnings for undesirable situations
185 
186  void PrintStatus(const G4double* StartArr,
187  G4double xstart,
188  const G4double* CurrentArr,
189  G4double xcurrent,
190  G4double requestStep,
191  G4int subStepNo);
192  void PrintStatus(const G4FieldTrack& StartFT,
193  const G4FieldTrack& CurrentFT,
194  G4double requestStep,
195  G4int subStepNo);
196  void PrintStat_Aux(const G4FieldTrack& aFieldTrack,
197  G4double requestStep,
198  G4double actualStep,
199  G4int subStepNo,
200  G4double subStepSize,
201  G4double dotVelocities);
202  // Verbose output for debugging
203 
204  void PrintStatisticsReport();
205  // Report on the number of steps, maximum errors etc.
206 
207 #ifdef QUICK_ADV_TWO
208  G4bool QuickAdvance( G4double yarrin[], // In
209  const G4double dydx[],
210  G4double hstep,
211  G4double yarrout[], // Out
212  G4double& dchord_step, // Out
213  G4double& dyerr ); // in length
214 #endif
215 
216  private:
217 
218  // ---------------------------------------------------------------
219  // INVARIANTS
220 
221  G4double fMinimumStep = 0.0;
222  // Minimum Step allowed in a Step (in absolute units)
223  G4double fSmallestFraction = 1.0e-12; // Expected range 1e-12 to 5e-15
224  // Smallest fraction of (existing) curve length - in relative units
225  // below this fraction the current step will be the last
226 
227  const G4int fNoIntegrationVariables = 0; // Variables in integration
228  const G4int fMinNoVars = 12; // Minimum number for FieldTrack
229  const G4int fNoVars = 0; // Full number of variable
230 
231  G4int fMaxNoSteps;
232  G4int fMaxStepBase = 250; // was 5000
233  // Default maximum number of steps is Base divided by the order of Stepper
234 
235  G4double safety;
236  G4double pshrnk; // exponent for shrinking
237  G4double pgrow; // exponent for growth
238  G4double errcon;
239  // Parameters used to grow and shrink trial stepsize.
240 
241  G4int fStatisticsVerboseLevel = 0;
242 
243  // ---------------------------------------------------------------
244  // DEPENDENT Objects
245 
246  G4MagIntegratorStepper* pIntStepper = nullptr;
247 
248  // ---------------------------------------------------------------
249  // STATE
250 
251  unsigned long fNoTotalSteps=0, fNoBadSteps=0;
252  unsigned long fNoSmallSteps=0, fNoInitialSmallSteps=0, fNoCalls=0;
253  G4double fDyerr_max=0.0, fDyerr_mx2=0.0;
254  G4double fDyerrPos_smTot=0.0, fDyerrPos_lgTot=0.0, fDyerrVel_lgTot=0.0;
255  G4double fSumH_sm=0.0, fSumH_lg=0.0;
256  // Step Statistics
257 
258  G4int fVerboseLevel = 0; // Verbosity level for printing (debug, ..)
259  // Could be varied during tracking - to help identify issues
260 
261  using ChordFinderDelegate = G4ChordFinderDelegate<G4MagInt_Driver>;
262 };
263 
264 #include "G4MagIntegratorDriver.icc"
265 
266 #endif