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