include/Geant4/G4BFieldIntegrationDriver.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 //
0027 //
0028 //
0029 // class G4BFieldIntegrationDriver
0030 //
0031 // Class description:
0032 //
0033 // Specialized integration driver for pure magnetic field
0034
0035 // History:
0036 // - Created. D.Sorokin
0037 // --------------------------------------------------------------------
0038
0039
0040 #include "globals.hh"
0041 #include "G4GeometryTolerance.hh"
0042 #include "G4FieldTrack.hh"
0043 #include "G4FieldUtils.hh"
0044
0045 namespace internal {
0046
0047 G4Mag_EqRhs* toMagneticEquation(G4EquationOfMotion* equation)
0048 {
0049     auto e = dynamic_cast<G4Mag_EqRhs*>(equation);
0050
0051     if (!e) {
0052         G4Exception("G4BFieldIntegrationDriver::G4BFieldIntegrationDriver",
0053                     "GeomField0003", FatalErrorInArgument,
0054                     "Works only with G4Mag_EqRhs");
0055     }
0056
0057     return e;
0058 }
0059
0060 } // internal
0061
0062
0063 template <class T>
0064 G4BFieldIntegrationDriver<T>::G4BFieldIntegrationDriver(G4double hminimum,
0065                                                         T* pStepper,
0066                                                         G4int numComponents,
0067                                                         G4int statisticsVerbose)
0068     : G4IntegrationDriver<T>(hminimum, pStepper, numComponents, statisticsVerbose)
0069     , fallbackThreshold(pi / 3.)
0070     , fequation(internal::toMagneticEquation(pStepper->GetEquationOfMotion()))
0071     , fallbackStepper(fequation)
0072 {
0073 }
0074
0075 template <class T>
0076 bool G4BFieldIntegrationDriver<T>::QuickAdvance(G4FieldTrack& fieldTrack,
0077                                                 const G4double dydx[],
0078                                                 G4double hstep,
0079                                                 G4double inverseCurvatureRadius,
0080                                                 G4double& dchord_step,
0081                                                 G4double& dyerr)
0082 {
0083     if (hstep * inverseCurvatureRadius < fallbackThreshold) {
0084         return G4IntegrationDriver<T>::QuickAdvance(
0085             fieldTrack, dydx, hstep, inverseCurvatureRadius, dchord_step, dyerr);
0086     }
0087
0088     G4IntegrationDriver<T>::IncrementQuickAdvanceCalls();
0089
0090     G4double yError[G4FieldTrack::ncompSVEC],
0091              yIn[G4FieldTrack::ncompSVEC],
0092              yOut[G4FieldTrack::ncompSVEC];
0093
0094     fieldTrack.DumpToArray(yIn);
0095
0096     fallbackStepper.Stepper(yIn, dydx, hstep, yOut, yError);
0097     dchord_step = fallbackStepper.DistChord();
0098     dyerr = field_utils::absoluteError(yOut, yError, hstep);
0099
0100     fieldTrack.LoadFromArray(yOut, fallbackStepper.GetNumberOfVariables());
0101     fieldTrack.SetCurveLength(fieldTrack.GetCurveLength() + hstep);
0102
0103     return true;
0104 }
0105
0106 template <class T>
0107 void G4BFieldIntegrationDriver<T>::SetEquationOfMotion(G4EquationOfMotion* equation)
0108 {
0109     G4IntegrationDriver<T>::SetEquationOfMotion(equation);
0110     fequation = internal::toMagneticEquation(equation);
0111 }
0112
0113 template <class T>
0114 G4double G4BFieldIntegrationDriver<T>::GetInverseCurvatureRadius(const G4FieldTrack& track,
0115                                                                  G4double field[]) const
0116 {
0117     const G4double Bmag = std::sqrt(field[0] * field[0] + field[1] * field[1] + field[2] * field[2]);
0118     const G4double momentum = track.GetMomentum().mag();
0119     const G4double particleCharge = fequation->FCof() / (CLHEP::eplus * CLHEP::c_light);
0120
0121     return std::abs(field_utils::inverseCurvatureRadius(particleCharge, momentum, Bmag));
0122 }