include/Geant4/G4QSStepper.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 // G4QSStepper
0027 //
0028 // QSS Integrator Stepper
0029 //
0030 // Authors: version 1 - Lucio Santi, Rodrigo Castro (Univ. Buenos Aires), 2018-2021
0031 //          version 2 - Mattias Portnoy (Univ. Buenos Aires), 2024
0032 // --------------------------------------------------------------------
0033 #ifndef G4QSS_STEPPER_HH
0034 #define G4QSS_STEPPER_HH
0035
0036 #include "G4FieldTrack.hh"
0037 #include "G4MagIntegratorStepper.hh"
0038 #include "G4QSSubstepStruct.hh"
0039
0040 #include <cmath>
0041 #include <CLHEP/Units/PhysicalConstants.h>
0042
0043 /**
0044  * @brief G4QSStepper is an integrator of particle's equation of
0045  * motion based on the QSS implementation.
0046  */
0047
0048 class G4QSStepper : public G4MagIntegratorStepper
0049 {
0050   public:
0051
0052     /**
0053      * Constructor for G4QSStepper.
0054      *  @param[in] equation Pointer to the provided equation of motion.
0055      *  @param[in] num_integration_vars The number of integration variables.
0056      *  @param[in] qssOrder The QSS order (2 or 3 expected; if <= 0 , use value
0057      *             from Messenger.
0058      */
0059     G4QSStepper( G4EquationOfMotion* equation,
0060                  G4int num_integration_vars = 6,  // always 6 -- ignore
0061                  G4int qssOrder= -1 );
0062
0063     /**
0064      * Default Destructor. Freeing of memory is done in susbsteps destructor.
0065      */
0066     ~G4QSStepper() override = default;
0067
0068     /**
0069      * Utility methods.
0070      */
0071     inline constexpr G4double Cubic_Function(const QSStateVector* states,
0072                                              G4int index, G4double delta_t);
0073     inline constexpr G4double Parabolic_Function(const QSStateVector* states,
0074                                                  G4int index, G4double delta_t);
0075     inline constexpr G4double Linear_Function(const QSStateVector* states,
0076                                               G4int index, G4double delta_t);
0077
0078     /**
0079      * 0 means position type, 1 means velocity type.
0080      */
0081     inline constexpr int INDEX_TYPE(G4int i);
0082
0083     /**
0084      * Auxiliary methods.
0085      */
0086     inline void momentum_to_velocity(const G4double* momentum, G4double* out);
0087     void set_relativistic_coeff(const G4double* momentum);
0088     inline void velocity_to_momentum(G4double *y);
0089
0090     /**
0091      * Key methods.
0092      */
0093     void initialize(const G4double y[]);
0094     inline void compare_time_and_update(G4int& index, G4int i);
0095     inline G4int get_next_sync_index();
0096     inline void update_field();
0097     inline G4double extrapolate_polynomial(QSStateVector* states,
0098                                     G4int index, G4double delta_t, G4int order);
0099     inline void extrapolate_all_states_to_t(Substep* substep,
0100                                             G4double t, G4double* yOut);
0101
0102     /**
0103      * Moves all the x states of variable index to the current time t.
0104      */
0105     inline void update_x(G4int index, G4double t);
0106
0107     /**
0108      * Moves all the q states of variable index to the current t.
0109      */
0110     inline void update_q(G4int index, G4double t);
0111
0112     /**
0113      * Update methods.
0114      */
0115     inline void update_x_position_derivates_using_q(G4int index);
0116     inline void update_x_velocity_derivates_using_q(G4int index);
0117     inline void update_x_derivates_using_q(G4int index);
0118     inline void update_sync_time_one_coefficient(G4int index);
0119
0120     /*
0121      * Updates when does the x,q distance goes beyond the quantum.
0122      * Uses polynomial roots-finding formulas.
0123      */
0124     void update_sync_time(G4int index);
0125
0126     /**
0127      * The stepper for the integration.
0128      * The stepsize is fixed, with the step size given by 'h'.
0129      * Integrates ODE starting values y[0 to 6]. Outputs yout[].
0130      *  @param[in] y Starting values array of integration variables.
0131      *  @param[in] dydx Derivatives array - Not used.
0132      *  @param[in] h The given step size.
0133      *  @param[out] yout Integration output.
0134      *  @param[out] yError The estimated error - Not used.
0135      */
0136     void Stepper( const G4double y[],
0137                   const G4double /*dydx*/ [],
0138                   G4double h,
0139                   G4double yout[],
0140                   G4double /* yerr */ [] ) override;
0141
0142     /**
0143      * Returns the QSS order of integration.
0144      */
0145     inline G4int IntegratorOrder() const override;
0146
0147     /**
0148      * Returns the stepper type-ID, "kQSStepper".
0149      */
0150     inline G4StepperType StepperType() const override { return kQSStepper; }
0151
0152     /**
0153      * Returns a pointer to the equation of motion.
0154      */
0155     inline G4EquationOfMotion* GetSpecificEquation();
0156
0157     /**
0158      * Returns current track state.
0159      */
0160     inline const field_utils::State& GetYOut() const;
0161
0162     /**
0163      * Track interpolation.
0164      *  @param[in] tau Step start, x.
0165      *  @param[in,out] yOut The current track state, y.
0166      */
0167     void Interpolate(G4double tau, G4double yOut[]);
0168
0169     /**
0170      * Returns the distance from chord line.
0171      */
0172     inline G4double DistChord() const override;
0173
0174     /**
0175      * Wrapper for the Stepper() function above.
0176      */
0177     inline void Stepper(const G4double yInput[],
0178                         const G4double dydx[],
0179                         G4double hstep,
0180                         G4double yOutput[],
0181                         G4double yError[],
0182                         G4double /*dydxOutput*/ []);
0183
0184     /**
0185      * Sets up interpolation. Does nothing.
0186      */
0187     inline void SetupInterpolation();
0188
0189     /*
0190      * Obligatory qss driver methods.
0191      */
0192     inline void reset(const G4FieldTrack* track);
0193     inline void SetPrecision(G4double dq_rel, G4double dq_min);
0194     inline G4double GetLastStepLength();
0195
0196     /*
0197      * Sets the mass at rest. Checking/ensuring that it is positive.
0198      */
0199     inline void setRestMass(G4double restMass);
0200
0201   private:
0202
0203     // Constants
0204
0205     static constexpr int DERIVATIVE_0{0};
0206     static constexpr int DERIVATIVE_1{1};
0207     static constexpr int DERIVATIVE_2{2};
0208     static constexpr int DERIVATIVE_3{3};
0209
0210     static constexpr int VX{3};
0211     static constexpr int VY{4};
0212     static constexpr int VZ{5};
0213
0214     static constexpr int POSITION_IDX{0};
0215     static constexpr int VELOCITY_IDX{3};
0216
0217     static constexpr G4double INFTY{1e+20};
0218
0219     /** Used to check if field changed from last update field during substeps. */
0220     G4bool fField_changed{true};
0221     G4bool fTrack_changed{true};
0222
0223     const G4int qss_order{2};
0224
0225     Substeps substeps;
0226     Substep current_substep;
0227     const G4FieldTrack* fCurrent_track{nullptr};
0228     QSStateVector dq_vector;
0229
0230     /** Invariants for this track -- during propagation. */
0231     G4double fCharge{-1.0};
0232     G4double fCharge_c2;
0233     G4double fRestMass{CLHEP::electron_mass_c2};
0234     G4double fGamma{1.0};
0235     G4double fCoeff; // coeff;
0236
0237     /** Cached values -- for tiny speed up. */
0238     G4double fMassOverC ; // was mass_times_gamma_over_speed_of_light;
0239     G4double fInv_mass_over_c;
0240
0241     /** Used by interpolation driver, need to copy state here when stepper finished. */
0242     G4double fYout[12];
0243
0244     /** QSS parameters separated into velocity and position. */
0245     G4double dqrel[2] = {0.0,0.0};
0246     G4double dqmin[2] = {0.001,0.001};
0247
0248     G4double fVelocity{0.0};
0249     G4double fFinal_t{0.0};
0250 };
0251
0252 // ----------------------------------------------------------------------------
0253 // Inline methods
0254 // ----------------------------------------------------------------------------
0255
0256 #include "G4QSStepper.icc"
0257
0258 #endif