include/Geant4/G4QSStepper.icc

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 inline methods implementation
0027 //
0028 // Authors: version 1 - Lucio Santi, Rodrigo Castro (Univ. Buenos Aires), 2018-2021
0029 //          version 2 - Mattias Portnoy (Univ. Buenos Aires), 2024
0030 // ----------------------------------------------------------------------------
0031
0032 inline
0033 constexpr G4double G4QSStepper::Cubic_Function(const QSStateVector* states,
0034                                                G4int index, G4double delta_t)
0035 {
0036   return states[DERIVATIVE_0][index] + (states[DERIVATIVE_1][index] + states[DERIVATIVE_2][index] * delta_t / 2 + states[DERIVATIVE_3][index] * delta_t * delta_t / 6) * delta_t;
0037 }
0038
0039 // ----------------------------------------------------------------------------
0040
0041 inline
0042 constexpr G4double G4QSStepper::Parabolic_Function(const QSStateVector* states,
0043                                                G4int index, G4double delta_t)
0044 {
0045   return states[DERIVATIVE_0][index] + (states[DERIVATIVE_1][index] + states[DERIVATIVE_2][index] * delta_t / 2) * delta_t;
0046 }
0047
0048 // ----------------------------------------------------------------------------
0049
0050 inline
0051 constexpr G4double G4QSStepper::Linear_Function(const QSStateVector* states,
0052                                                G4int index, G4double delta_t)
0053 {
0054   return states[DERIVATIVE_0][index] + states[DERIVATIVE_1][index] * delta_t;
0055 }
0056
0057 // ----------------------------------------------------------------------------
0058
0059 inline
0060 constexpr G4int G4QSStepper::INDEX_TYPE(G4int i)
0061 {
0062   return i >> 2;
0063 }
0064
0065 // ----------------------------------------------------------------------------
0066
0067 inline
0068 void G4QSStepper::momentum_to_velocity(const G4double* momentum, G4double* out)
0069 {
0070   out[0] = momentum[0] * fInv_mass_over_c;
0071   out[1] = momentum[1] * fInv_mass_over_c;
0072   out[2] = momentum[2] * fInv_mass_over_c;
0073 }
0074
0075 // ----------------------------------------------------------------------------
0076
0077 inline
0078 void G4QSStepper::compare_time_and_update(G4int& index, G4int i)
0079 {
0080   if (current_substep.sync_t[i] < current_substep.sync_t[index]) { index = i;}
0081 }
0082
0083 // ----------------------------------------------------------------------------
0084
0085 inline
0086 G4int G4QSStepper::IntegratorOrder() const
0087 {
0088   return qss_order;
0089 }
0090
0091 // ----------------------------------------------------------------------------
0092
0093 inline
0094 G4EquationOfMotion* G4QSStepper::GetSpecificEquation()
0095 {
0096   return GetEquationOfMotion();
0097 }
0098
0099 // ----------------------------------------------------------------------------
0100
0101 inline
0102 const field_utils::State& G4QSStepper::GetYOut() const
0103 {
0104   return fYout;
0105 }
0106
0107 // ----------------------------------------------------------------------------
0108
0109 inline
0110 G4double G4QSStepper::DistChord() const
0111 {
0112   return 0.;
0113 }
0114
0115 // ----------------------------------------------------------------------------
0116
0117 inline
0118 void G4QSStepper::Stepper(const G4double yInput[],
0119              const G4double dydx[], G4double hstep,
0120              G4double yOutput[], G4double yError[], G4double /*dydxOutput*/ [])
0121 {
0122   Stepper(yInput, dydx, hstep, yOutput, yError);
0123 }
0124
0125 // ----------------------------------------------------------------------------
0126
0127 inline
0128 void G4QSStepper::SetupInterpolation()
0129 {
0130 }
0131
0132 inline
0133 void G4QSStepper::reset(const G4FieldTrack *track)
0134 {
0135   fTrack_changed = true; //// Cannot rely on addresses --- OLD was track != fCurrent_track;
0136   fCurrent_track = track;
0137   setRestMass(track->GetRestMass());
0138 }
0139
0140 // ----------------------------------------------------------------------------
0141
0142 inline
0143 void G4QSStepper::setRestMass(G4double restMass)
0144 {
0145   if( restMass > 0.0 ){
0146       fRestMass = restMass;
0147   } else {
0148       fRestMass = 1.0e-10 * CLHEP::electronvolt;
0149 #ifdef G4DEBUG_FIELD
0150       G4cerr << "G4QSStepper::setRestMass(): WARNING: Rest Mass (track) = " << restMass << " Reset to:" << fRestMass / CLHEP::eV << " eV " << G4endl;
0151 #endif
0152   }
0153 }
0154
0155 // ----------------------------------------------------------------------------
0156
0157 inline
0158 void G4QSStepper::SetPrecision(G4double dq_rel, G4double dq_min)
0159 {
0160   dqmin[0] = dq_min;
0161   dqmin[1] = dq_min;
0162   dqrel[0] = dq_rel;
0163   dqrel[1] = dq_rel;
0164 }
0165
0166 // ----------------------------------------------------------------------------
0167
0168 inline
0169 G4double G4QSStepper::GetLastStepLength()
0170 {
0171   return current_substep.t * fVelocity;
0172 }
0173
0174 // ----------------------------------------------------------------------------
0175
0176 inline
0177 void G4QSStepper::velocity_to_momentum(G4double *y)
0178 {
0179   y[3] *= fMassOverC;
0180   y[4] *= fMassOverC;
0181   y[5] *= fMassOverC;
0182 }
0183
0184 // ----------------------------------------------------------------------------
0185
0186 inline
0187 G4int G4QSStepper::get_next_sync_index()
0188 {
0189
0190   // Goes through each index and get the one with the closest sync t.
0191   // Unrolled loop for tiny speedup.
0192
0193   G4int index = 0;
0194   compare_time_and_update(index,1);
0195   compare_time_and_update(index,2);
0196   compare_time_and_update(index,3);
0197   compare_time_and_update(index,4);
0198   compare_time_and_update(index,5);
0199   return index;
0200 }
0201
0202 // ----------------------------------------------------------------------------
0203
0204 inline
0205 void G4QSStepper::update_field()
0206 {
0207   G4double old_field[3] = { current_substep.b_field[0],
0208                             current_substep.b_field[1],
0209                             current_substep.b_field[2] };
0210   GetEquationOfMotion()->GetFieldValue(current_substep.state_x[DERIVATIVE_0],
0211                                        current_substep.b_field );
0212   fField_changed = false;
0213   for (G4int i = 0; i < 3 && ! fField_changed; ++i)
0214   {
0215     fField_changed = fField_changed || old_field[i] != current_substep.b_field[i];
0216   }
0217 }
0218
0219 // ----------------------------------------------------------------------------
0220
0221 inline
0222 G4double G4QSStepper::extrapolate_polynomial(QSStateVector* states,
0223                                 G4int index, G4double delta_t, G4int order)
0224 {
0225   if (delta_t == 0 || order == 0) { return states[DERIVATIVE_0][index]; }
0226
0227   switch (order)
0228   {
0229     case 2:
0230       return Parabolic_Function(states,index,delta_t);
0231       break;
0232     case 3:
0233       return Cubic_Function(states,index,delta_t);
0234       break;
0235     case 1:
0236       return Linear_Function(states,index,delta_t);
0237       break;
0238     default:
0239       // TODO check how to raise error
0240       return 146546;
0241   }
0242 }
0243
0244 // ----------------------------------------------------------------------------
0245
0246 inline
0247 void G4QSStepper::extrapolate_all_states_to_t(Substep* substep,
0248                                               G4double t, G4double* yOut)
0249 {
0250   for (G4int j = 0; j < 6; ++j)
0251   {
0252     G4double t_j = substep->state_tx[j];
0253     G4double delta_tj = t - t_j;
0254     yOut[j] = extrapolate_polynomial(&substep->state_x[DERIVATIVE_0], j, delta_tj, substep->extrapolation_method);
0255   }
0256 }
0257
0258 // ----------------------------------------------------------------------------
0259
0260 inline
0261 void G4QSStepper::update_x(G4int index, G4double t)
0262 {
0263   G4double delta_t = t - current_substep.state_tx[index];
0264   if (delta_t == 0) { return; }
0265
0266   //current_substep.state_x[DERIVATE_1][index] += current_substep.state_x[DERIVATE_2][index] * delta_t;
0267   switch (qss_order)
0268   {
0269     case 2:
0270       current_substep.state_x[DERIVATIVE_0][index] = Parabolic_Function(current_substep.state_x,index,delta_t);
0271       current_substep.state_x[DERIVATIVE_1][index] = Linear_Function((&current_substep.state_x[DERIVATIVE_1]),index,delta_t);
0272       break;
0273     case 3:
0274       current_substep.state_x[DERIVATIVE_0][index] = Cubic_Function(current_substep.state_x,index,delta_t);
0275       current_substep.state_x[DERIVATIVE_1][index] = Parabolic_Function((&current_substep.state_x[DERIVATIVE_1]),index,delta_t);
0276       current_substep.state_x[DERIVATIVE_2][index] = Linear_Function((&current_substep.state_x[DERIVATIVE_2]),index,delta_t);
0277       break;
0278     case 1:
0279       current_substep.state_x[DERIVATIVE_0][index] = Linear_Function(current_substep.state_x,index,delta_t);
0280       break;
0281     default:
0282       break;
0283   }
0284   current_substep.state_tx[index] = t;
0285 }
0286
0287 // ----------------------------------------------------------------------------
0288
0289 inline
0290 void G4QSStepper::update_q(G4int index, G4double t)
0291 {
0292   G4double delta_t = t - current_substep.state_tq[index];
0293   if (delta_t == 0) { return; }
0294   switch (qss_order)
0295   {
0296     case 2:
0297       current_substep.state_q[DERIVATIVE_0][index] = Linear_Function(current_substep.state_q,index,delta_t);
0298       break;
0299     case 3:
0300       current_substep.state_q[DERIVATIVE_0][index] = Parabolic_Function(current_substep.state_q,index,delta_t);
0301       current_substep.state_q[DERIVATIVE_1][index] = Linear_Function((&current_substep.state_q[DERIVATIVE_1]),index,delta_t);
0302       break;
0303     case 1:
0304       break;
0305   }
0306   current_substep.state_tq[index] = t;
0307
0308 }
0309
0310 // ----------------------------------------------------------------------------
0311
0312 inline
0313 void G4QSStepper::update_x_position_derivates_using_q(G4int index)
0314 {
0315   // assumes index is position index
0316
0317   current_substep.state_x[DERIVATIVE_1][index] =
0318       current_substep.state_q[DERIVATIVE_0][index+VELOCITY_IDX];
0319   current_substep.state_x[DERIVATIVE_2][index] =
0320       current_substep.state_q[DERIVATIVE_1][index+VELOCITY_IDX];
0321   current_substep.state_x[DERIVATIVE_3][index] =
0322       current_substep.state_q[DERIVATIVE_2][index+VELOCITY_IDX];
0323 }
0324
0325 // ----------------------------------------------------------------------------
0326
0327 inline
0328 void G4QSStepper::update_x_velocity_derivates_using_q(G4int index)
0329 {
0330   // asumes index is velocity index
0331
0332   G4int modulo = VELOCITY_IDX;
0333   G4int index_pos = (index+modulo+1)%modulo;
0334   G4int index_neg = (index+modulo-1)%modulo;
0335
0336   G4double b1 = current_substep.b_field[index_pos];
0337   G4double b2 = current_substep.b_field[index_neg];
0338   for (G4int derivate_order = 0; derivate_order < qss_order; ++derivate_order)
0339   {
0340     current_substep.state_x[derivate_order+1][index] =
0341       fCoeff* (
0342       current_substep.state_q[derivate_order][index_pos+VELOCITY_IDX] * b2 -
0343       current_substep.state_q[derivate_order][index_neg+VELOCITY_IDX] * b1
0344       );
0345   }
0346 }
0347
0348 // ----------------------------------------------------------------------------
0349
0350 inline
0351 void G4QSStepper::update_x_derivates_using_q(G4int index)
0352 {
0353   // updates x using q with the Lorentz equation
0354
0355   if (index < VELOCITY_IDX)
0356   {
0357     update_x_position_derivates_using_q(index);
0358   }
0359   else
0360   {
0361     update_x_velocity_derivates_using_q(index);
0362   }
0363 }
0364
0365 // ----------------------------------------------------------------------------
0366
0367 /*  Updates when does the x,q distance goes beyond the quantum.
0368     For the special case of both polynomials being equal except
0369     for higher coefficient- Such as after syncing */
0370 inline
0371 void G4QSStepper::update_sync_time_one_coefficient(G4int index)
0372 {
0373   G4double leading_poly_cofficient = current_substep.state_x[qss_order][index];
0374   if (leading_poly_cofficient == 0)
0375   {
0376       current_substep.sync_t[index] = INFTY;
0377   }
0378   else
0379   {
0380     G4double dq_leading_ratio = dq_vector[index]/std::fabs(leading_poly_cofficient);
0381     switch (qss_order)
0382     {
0383       case 2:
0384         current_substep.sync_t[index] = current_substep.state_tx[index] + std::sqrt(dq_leading_ratio);
0385         break;
0386       case 3:
0387         current_substep.sync_t[index] = current_substep.state_tx[index] + cbrt(dq_leading_ratio);
0388         break;
0389       case 1:
0390         current_substep.sync_t[index] = current_substep.state_tx[index] + dq_leading_ratio;
0391         break;
0392     }
0393   }
0394 }