Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4NystromRK4.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// G4NystromRK4 implmentation
27//
28// Created: I.Gavrilenko, 15.05.2009 (as G4AtlasRK4)
29// Adaptations: J.Apostolakis, November 2009
30// -------------------------------------------------------------------
31
32#include <memory>
33
34#include "G4NystromRK4.hh"
35
36#include "G4Exception.hh"
37#include "G4SystemOfUnits.hh"
38#include "G4FieldUtils.hh"
39#include "G4LineSection.hh"
40
41using namespace field_utils;
42
43namespace
44{
45 G4bool notEquals(G4double p1, G4double p2)
46 {
47 return std::fabs(p1 - p2) > perMillion * p2;
48 }
49 constexpr G4int INTEGRATED_COMPONENTS = 6;
50} // namespace
51
52
53G4NystromRK4::G4NystromRK4(G4Mag_EqRhs* equation, G4double distanceConstField)
54 : G4MagIntegratorStepper(equation, INTEGRATED_COMPONENTS)
55{
56 if (distanceConstField > 0)
57 {
58 SetDistanceForConstantField(distanceConstField);
59 }
60}
61
63 const G4double dydx[],
64 G4double hstep,
65 G4double yOut[],
66 G4double yError[])
67{
68 fInitialPoint = { y[0], y[1], y[2] };
69
70 G4double field[3];
71
72 constexpr G4double one_sixth= 1./6.;
73 const G4double S5 = 0.5 * hstep;
74 const G4double S4 = .25 * hstep;
75 const G4double S6 = hstep * one_sixth;
76
77 const G4double momentum2 = getValue2(y, Value3D::Momentum);
78 if (notEquals(momentum2, fMomentum2))
79 {
80 fMomentum = std::sqrt(momentum2);
81 fMomentum2 = momentum2;
82 fInverseMomentum = 1. / fMomentum;
83 fCoefficient = GetFCof() * fInverseMomentum;
84 }
85
86 // Point 1
87 const G4double K1[3] = {
88 fInverseMomentum * dydx[3],
89 fInverseMomentum * dydx[4],
90 fInverseMomentum * dydx[5]
91 };
92
93 // Point2
94 G4double p[4] = {
95 y[0] + S5 * (dydx[0] + S4 * K1[0]),
96 y[1] + S5 * (dydx[1] + S4 * K1[1]),
97 y[2] + S5 * (dydx[2] + S4 * K1[2]),
98 y[7]
99 };
100
101 GetFieldValue(p, field);
102
103 const G4double A2[3] = {
104 dydx[0] + S5 * K1[0],
105 dydx[1] + S5 * K1[1],
106 dydx[2] + S5 * K1[2]
107 };
108
109 const G4double K2[3] = {
110 (A2[1] * field[2] - A2[2] * field[1]) * fCoefficient,
111 (A2[2] * field[0] - A2[0] * field[2]) * fCoefficient,
112 (A2[0] * field[1] - A2[1] * field[0]) * fCoefficient
113 };
114
115 fMidPoint = { p[0], p[1], p[2] };
116
117 // Point 3 with the same magnetic field
118 const G4double A3[3] = {
119 dydx[0] + S5 * K2[0],
120 dydx[1] + S5 * K2[1],
121 dydx[2] + S5 * K2[2]
122 };
123
124 const G4double K3[3] = {
125 (A3[1] * field[2] - A3[2] * field[1]) * fCoefficient,
126 (A3[2] * field[0] - A3[0] * field[2]) * fCoefficient,
127 (A3[0] * field[1] - A3[1] * field[0]) * fCoefficient
128 };
129
130 // Point 4
131 p[0] = y[0] + hstep * (dydx[0] + S5 * K3[0]);
132 p[1] = y[1] + hstep * (dydx[1] + S5 * K3[1]);
133 p[2] = y[2] + hstep * (dydx[2] + S5 * K3[2]);
134
135 GetFieldValue(p, field);
136
137 const G4double A4[3] = {
138 dydx[0] + hstep * K3[0],
139 dydx[1] + hstep * K3[1],
140 dydx[2] + hstep * K3[2]
141 };
142
143 const G4double K4[3] = {
144 (A4[1] * field[2] - A4[2] * field[1]) * fCoefficient,
145 (A4[2] * field[0] - A4[0] * field[2]) * fCoefficient,
146 (A4[0] * field[1] - A4[1] * field[0]) * fCoefficient
147 };
148
149 // New position
150 yOut[0] = y[0] + hstep * (dydx[0] + S6 * (K1[0] + K2[0] + K3[0]));
151 yOut[1] = y[1] + hstep * (dydx[1] + S6 * (K1[1] + K2[1] + K3[1]));
152 yOut[2] = y[2] + hstep * (dydx[2] + S6 * (K1[2] + K2[2] + K3[2]));
153 // New direction
154 yOut[3] = dydx[0] + S6 * (K1[0] + K4[0] + 2. * (K2[0] + K3[0]));
155 yOut[4] = dydx[1] + S6 * (K1[1] + K4[1] + 2. * (K2[1] + K3[1]));
156 yOut[5] = dydx[2] + S6 * (K1[2] + K4[2] + 2. * (K2[2] + K3[2]));
157 // Pass Energy, time unchanged -- time is not integrated !!
158 yOut[6] = y[6];
159 yOut[7] = y[7];
160
161 fEndPoint = { yOut[0], yOut[1], yOut[2] };
162
163 // Errors
164 yError[3] = hstep * std::fabs(K1[0] - K2[0] - K3[0] + K4[0]);
165 yError[4] = hstep * std::fabs(K1[1] - K2[1] - K3[1] + K4[1]);
166 yError[5] = hstep * std::fabs(K1[2] - K2[2] - K3[2] + K4[2]);
167 yError[0] = hstep * yError[3];
168 yError[1] = hstep * yError[4];
169 yError[2] = hstep * yError[5];
170 yError[3] *= fMomentum;
171 yError[4] *= fMomentum;
172 yError[5] *= fMomentum;
173
174 // Normalize momentum
175 const G4double normF = fMomentum / getValue(yOut, Value3D::Momentum);
176 yOut[3] *= normF;
177 yOut[4] *= normF;
178 yOut[5] *= normF;
179
180 // My trial code:
181 // G4double endMom2 = yOut[3]*yOut[3]+yOut[4]*yOut[4]+yOut[5]*yOut[5];
182 // G4double normF = std::sqrt( startMom2 / endMom2 );
183}
184
186{
187 return G4LineSection::Distline(fMidPoint, fInitialPoint, fEndPoint);
188}
189
191{
192 if (GetField() == nullptr)
193 {
194 G4Exception("G4NystromRK4::SetDistanceForConstantField",
195 "Nystrom 001", JustWarning,
196 "Provided field is not G4CachedMagneticField. Changing field type.");
197
198 fCachedField = std::make_unique<G4CachedMagneticField>(
199 dynamic_cast<G4MagneticField*>(GetEquationOfMotion()->GetFieldObj()),
200 length);
201
202 GetEquationOfMotion()->SetFieldObj(fCachedField.get());
203 }
204 GetField()->SetConstDistance(length);
205}
206
208{
209 if (GetField() == nullptr)
210 {
211 return 0.0;
212 }
213 return GetField()->GetConstDistance();
214}
215
216G4CachedMagneticField* G4NystromRK4::GetField()
217{
218 return dynamic_cast<G4CachedMagneticField*>(GetEquationOfMotion()->GetFieldObj());
219}
220
221const G4CachedMagneticField* G4NystromRK4::GetField() const
222{
223 return const_cast<G4NystromRK4*>(this)->GetField();
224}
@ JustWarning
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
double G4double
Definition G4Types.hh:83
bool G4bool
Definition G4Types.hh:86
int G4int
Definition G4Types.hh:85
G4double GetConstDistance() const
void SetConstDistance(G4double dist)
const G4Field * GetFieldObj() const
void SetFieldObj(G4Field *pField)
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
G4EquationOfMotion * GetEquationOfMotion()
G4double GetDistanceForConstantField() const
void Stepper(const G4double y[], const G4double dydx[], G4double hstep, G4double yOut[], G4double yError[]) override
G4double DistChord() const override
void SetDistanceForConstantField(G4double length)
G4NystromRK4(G4Mag_EqRhs *EquationMotion, G4double distanceConstField=0.0)
G4double getValue(const ArrayType &array, Value1D value)
G4double getValue2(const ArrayType &array, Value1D value)