include/Geant4/G4OldMagIntDriver.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // G4OldMagIntDriver
0027 //
0028 // Class description:
0029 //
0030 // Provides a driver that talks to the Integrator Stepper, and insures that
0031 // the error is within acceptable bounds.
0032
0033 // Author: Vladimir Grichine (CERN), 07.10.1996 - Created
0034 //         W.Wander (MIT), 28.01.1998 - Added ability for low order integrators
0035 // --------------------------------------------------------------------
0036 #ifndef G4OLD_MAGINT_DRIVER_HH
0037 #define G4OLD_MAGINT_DRIVER_HH
0038
0039 #include "G4VIntegrationDriver.hh"
0040 #include "G4MagIntegratorStepper.hh"
0041 #include "G4ChordFinderDelegate.hh"
0042
0043 /**
0044  * @brief G4OldMagIntDriver provides a driver that talks to the Integrator
0045  * Stepper and insures that the error is within acceptable bounds.
0046  */
0047
0048 class G4OldMagIntDriver : public G4VIntegrationDriver,
0049                           public G4ChordFinderDelegate<G4OldMagIntDriver>
0050 {
0051   public:
0052
0053     /**
0054      * Constructor for G4OldMagIntDriver.
0055      *  @param[in] hminimum The minumum allowed step.
0056      *  @param[in] pItsStepper Pointer to the integrator stepper.
0057      *  @param[in] numberOfComponents The number of integration variables.
0058      *  @param[in] statisticsVerbosity Flag for verbosity.
0059      */
0060     G4OldMagIntDriver(G4double hminimum,
0061                       G4MagIntegratorStepper* pItsStepper,
0062                       G4int numberOfComponents = 6,
0063                       G4int statisticsVerbosity = 0);
0064
0065     /**
0066      * Destructor. Provides statistics if verbosity level is greater than 1.
0067      */
0068     ~G4OldMagIntDriver() override;
0069
0070     /**
0071      * Copy constructor and assignment operator not allowed.
0072      */
0073     G4OldMagIntDriver(const G4OldMagIntDriver&) = delete;
0074     G4OldMagIntDriver& operator=(const G4OldMagIntDriver&) = delete;
0075
0076     /**
0077      * Computes the step to take, based on chord limits.
0078      *  @param[in,out] track The current track in field.
0079      *  @param[in] stepMax Proposed maximum step length.
0080      *  @param[in] epsStep Requested accuracy, y_err/hstep.
0081      *  @param[in] chordDistance Maximum sagitta distance.
0082      *  @returns The length of step taken.
0083      */
0084     inline G4double AdvanceChordLimited(G4FieldTrack& track,
0085                                         G4double stepMax,
0086                                         G4double epsStep,
0087                                         G4double chordDistance) override;
0088
0089     /**
0090      * Dispatch interface method for initialisation/reset of driver.
0091      */
0092     inline void OnStartTracking() override;
0093
0094     /**
0095      * Dispatch interface method for computing step. Does nothing here.
0096      */
0097     inline void OnComputeStep(const G4FieldTrack* = nullptr) override;
0098
0099     /**
0100      * The driver implements re-integration, so returns true.
0101      */
0102     inline G4bool DoesReIntegrate() const override;
0103
0104     /**
0105      * Advances integration accurately by relative accuracy better than 'eps'.
0106      *  @param[in,out] y_current The current track in field.
0107      *  @param[in] hstep Proposed step length.
0108      *  @param[in] eps Requested accuracy, y_err/hstep.
0109      *  @param[in] hinitial Initial minimum integration step.
0110      *  @returns true if integration succeeds.
0111      */
0112     G4bool AccurateAdvance(G4FieldTrack& y_current,
0113                            G4double hstep,
0114                            G4double eps,  // Requested y_err/hstep
0115                            G4double hinitial = 0.0) override;
0116
0117     /**
0118      * Attempts one integration step, and returns estimated error 'dyerr'.
0119      * It does not ensure accuracy.
0120      *  @param[in,out] y_val The current track in field.
0121      *  @param[in] dydx dydx array.
0122      *  @param[in] hstep Proposed step length.
0123      *  @param[out] dchord_step Estimated sagitta distance.
0124      *  @param[out] dyerr Estimated error.
0125      *  @returns true if integration succeeds.
0126      */
0127     G4bool QuickAdvance(G4FieldTrack& y_val,      // In/Out
0128                                 const G4double dydx[],
0129                                 G4double hstep,
0130                                 G4double& dchord_step,
0131                                 G4double& dyerr) override;
0132
0133     /**
0134      * Attempts one integration step, and returns estimated error 'dyerr'.
0135      * It does not ensure accuracy.
0136      *  @param[in,out] y_posvel The current track in field.
0137      *  @param[in] dydx dydx array.
0138      *  @param[in] hstep Proposed step length.
0139      *  @param[out] dchord_step Estimated sagitta distance.
0140      *  @param[out] dyerr_pos_sq Estimated error in position.
0141      *  @param[out] dyerr_mom_rel_sq Estimated error in momentum
0142      *              (normalised: Delta_Integration(p^2)/(p^2)).
0143      *  @returns true if integration succeeds.
0144      */
0145     G4bool QuickAdvance(G4FieldTrack& y_posvel,   // In/Out
0146                         const G4double dydx[],
0147                         G4double hstep,           // In
0148                         G4double& dchord_step,
0149                         G4double& dyerr_pos_sq,
0150                         G4double& dyerr_mom_rel_sq);
0151
0152     /**
0153      * Accessors.
0154      */
0155     inline G4double GetHmin() const;
0156     inline G4double Hmin() const;     // Obsolete
0157     inline G4double GetSafety() const;
0158     inline G4double GetPshrnk() const;
0159     inline G4double GetPgrow() const;
0160     inline G4double GetErrcon() const;
0161     void GetDerivatives(const G4FieldTrack& y_curr,            // INput
0162                               G4double dydx[]) const override; // OUTput
0163     void GetDerivatives(const G4FieldTrack& track,
0164                               G4double dydx[],
0165                               G4double field[]) const override;
0166     /**
0167      * Getter and setter for the equation of motion.
0168      */
0169     G4EquationOfMotion* GetEquationOfMotion() override;
0170     void SetEquationOfMotion(G4EquationOfMotion* equation) override;
0171
0172     /**
0173      * Sets a new stepper 'pItsStepper' for this driver. Then it calls
0174      * ResetParameters() to update its parameters accordingly.
0175      */
0176     void RenewStepperAndAdjust(G4MagIntegratorStepper* pItsStepper) override;
0177
0178     /**
0179      * Resets the qarameters according to the new provided safety value.
0180      *  i) sets the exponents (pgrow & pshrnk), using the current order;
0181      * ii) sets the safety and calculates "errcon" according to the above values.
0182      */
0183     inline void ReSetParameters(G4double new_safety = 0.9);
0184
0185     /**
0186      * Modifiers. When setting safety or pgrow, errcon will be set
0187      * to a compatible value.
0188      */
0189     inline void SetSafety(G4double valS);
0190     inline void SetPshrnk(G4double valPs);
0191     inline void SetPgrow (G4double valPg);
0192     inline void SetErrcon(G4double valEc);
0193     inline G4double ComputeAndSetErrcon();
0194
0195     /**
0196      * Accessors for the integrator stepper.
0197      */
0198     const G4MagIntegratorStepper* GetStepper() const override;
0199     G4MagIntegratorStepper* GetStepper() override;
0200
0201     /**
0202      * Takes one Step that is as large as possible while satisfying the
0203      * accuracy criterion of: yerr < eps * |y_end-y_start|.
0204      *  @param[in,out] ystart The current track state, y.
0205      *  @param[in] dydx The derivatives array.
0206      *  @param[in,out] x Step start, x.
0207      *  @param[in] htry Step to attempt.
0208      *  @param[in] eps The relative accuracy.
0209      *  @param[out] hdid Step achieved.
0210      *  @param[out] hnext Proposed next step.
0211      *  @returns true if integration succeeds.
0212      */
0213     void OneGoodStep(G4double ystart[], // Like old RKF45step()
0214                      const G4double dydx[],
0215                      G4double& x,
0216                      G4double htry,
0217                      G4double eps,
0218                      G4double& hdid,
0219                      G4double& hnext) ;
0220
0221     G4double ComputeNewStepSize(G4double errMaxNorm, // normalised
0222                                 G4double hstepCurrent) override;
0223       // Taking the last step's normalised error, calculate
0224       // a step size for the next step.
0225       // Do not limit the next step's size within a factor of the
0226       // current one.
0227
0228     /**
0229      * Writes out to stream the parameters/state of the driver.
0230      */
0231     void StreamInfo( std::ostream& os ) const override;
0232
0233     /**
0234      * Takes the last step's normalised error and calculates a step size
0235      * for the next step. Limits the next step's size within a range around
0236      * the current one.
0237      */
0238     G4double ComputeNewStepSize_WithinLimits(G4double errMaxNorm, // normalised
0239                                              G4double hstepCurrent);
0240
0241     /**
0242      * Modifier and accessor for the maximum number of steps that can be taken
0243      * for the integration of a single segment, i.e. a single call to
0244      * AccurateAdvance().
0245      */
0246     inline G4int GetMaxNoSteps() const;
0247     inline void SetMaxNoSteps(G4int val);
0248
0249     /**
0250      * More modifiers and accessors.
0251      */
0252     inline void SetHmin(G4double newval);
0253     void SetVerboseLevel(G4int newLevel) override;
0254     G4int GetVerboseLevel() const override;
0255     inline G4double GetSmallestFraction() const;
0256     void SetSmallestFraction( G4double val );
0257
0258   protected:
0259
0260     /**
0261      * Loggers, issuing warnings for undesirable situations.
0262      */
0263     void WarnSmallStepSize(G4double hnext, G4double hstep,
0264                            G4double h, G4double xDone,
0265                            G4int noSteps);
0266     void WarnTooManyStep(G4double x1start, G4double x2end, G4double xCurrent);
0267     void WarnEndPointTooFar(G4double endPointDist,
0268                             G4double hStepSize ,
0269                             G4double epsilonRelative,
0270                             G4int debugFlag);
0271
0272     /**
0273      * Loggers for verbosity printouts.
0274      */
0275     void PrintStatus(const G4double* StartArr,
0276                            G4double xstart,
0277                      const G4double* CurrentArr,
0278                            G4double xcurrent,
0279                            G4double requestStep,
0280                            G4int subStepNo);
0281     void PrintStatus(const G4FieldTrack& StartFT,
0282                      const G4FieldTrack& CurrentFT,
0283                            G4double requestStep,
0284                            G4int subStepNo);
0285     void PrintStat_Aux(const G4FieldTrack& aFieldTrack,
0286                              G4double requestStep,
0287                              G4double actualStep,
0288                              G4int subStepNo,
0289                              G4double subStepSize,
0290                              G4double dotVelocities);
0291
0292     /**
0293      * Reports on the number of steps, maximum errors etc.
0294      */
0295     void PrintStatisticsReport();
0296
0297 #ifdef QUICK_ADV_TWO
0298      G4bool QuickAdvance(      G4double  yarrin[],     // In
0299                          const G4double  dydx[],
0300                                G4double  hstep,
0301                                G4double  yarrout[],    // Out
0302                                G4double& dchord_step,  // Out
0303                                G4double& dyerr );      // in length
0304 #endif
0305
0306   private:
0307
0308     // ---------------------------------------------------------------
0309     //  INVARIANTS
0310
0311     /** Minimum Step allowed in a Step (in absolute units). */
0312     G4double fMinimumStep = 0.0;
0313
0314     /** Smallest fraction of (existing) curve length, in relative units.
0315         Below this fraction the current step will be the last. */
0316     G4double fSmallestFraction = 1.0e-12;    // Expected range 1e-12 to 5e-15
0317
0318     /** Variables in integration. */
0319     const G4int fNoIntegrationVariables = 0;
0320
0321     /** Minimum number for FieldTrack. */
0322     const G4int fMinNoVars = 12;
0323
0324     /** Full number of variable. */
0325     const G4int fNoVars = 0;
0326
0327     /** Default maximum number of steps is Base divided by the order of Stepper. */
0328     G4int fMaxNoSteps;
0329     G4int fMaxStepBase = 250;  // was 5000
0330
0331     /** Parameters used to grow and shrink trial stepsize. */
0332     G4double safety;
0333     G4double pshrnk;   //  exponent for shrinking
0334     G4double pgrow;    //  exponent for growth
0335     G4double errcon;
0336
0337     G4int fStatisticsVerboseLevel = 0;
0338
0339     // ---------------------------------------------------------------
0340     // DEPENDENT Objects
0341
0342     G4MagIntegratorStepper* pIntStepper = nullptr;
0343
0344     // ---------------------------------------------------------------
0345     //  STATE
0346
0347     /** Step Statistics. */
0348     unsigned long fNoTotalSteps=0, fNoBadSteps=0;
0349     unsigned long fNoSmallSteps=0, fNoInitialSmallSteps=0, fNoCalls=0;
0350     G4double fDyerr_max=0.0, fDyerr_mx2=0.0;
0351     G4double fDyerrPos_smTot=0.0, fDyerrPos_lgTot=0.0, fDyerrVel_lgTot=0.0;
0352     G4double fSumH_sm=0.0, fSumH_lg=0.0;
0353
0354     /** Could be varied during tracking - to help identify issues. */
0355     G4int fVerboseLevel = 0;   // Verbosity level for printing (debug, ..)
0356
0357     using ChordFinderDelegate = G4ChordFinderDelegate<G4OldMagIntDriver>;
0358 };
0359
0360 #include "G4OldMagIntDriver.icc"
0361
0362 #endif