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