Geant4 Cross Reference |
1 // ******************************************* 1 // ******************************************************************** 2 // * License and Disclaimer 2 // * License and Disclaimer * 3 // * 3 // * * 4 // * The Geant4 software is copyright of th 4 // * The Geant4 software is copyright of the Copyright Holders of * 5 // * the Geant4 Collaboration. It is provided 5 // * the Geant4 Collaboration. It is provided under the terms and * 6 // * conditions of the Geant4 Software License 6 // * conditions of the Geant4 Software License, included in the file * 7 // * LICENSE and available at http://cern.ch/ 7 // * LICENSE and available at http://cern.ch/geant4/license . These * 8 // * include a list of copyright holders. 8 // * include a list of copyright holders. * 9 // * 9 // * * 10 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 11 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 12 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 13 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 14 // * use. Please see the license in the file 14 // * use. Please see the license in the file LICENSE and URL above * 15 // * for the full disclaimer and the limitatio 15 // * for the full disclaimer and the limitation of liability. * 16 // * 16 // * * 17 // * This code implementation is the result 17 // * This code implementation is the result of the scientific and * 18 // * technical work of the GEANT4 collaboratio 18 // * technical work of the GEANT4 collaboration. * 19 // * By using, copying, modifying or distri 19 // * By using, copying, modifying or distributing the software (or * 20 // * any work based on the software) you ag 20 // * any work based on the software) you agree to acknowledge its * 21 // * use in resulting scientific publicati 21 // * use in resulting scientific publications, and indicate your * 22 // * acceptance of all terms of the Geant4 Sof 22 // * acceptance of all terms of the Geant4 Software license. * 23 // ******************************************* 23 // ******************************************************************** 24 // 24 // 25 // Helper namespace field_utils implementation 25 // Helper namespace field_utils implementation 26 // 26 // 27 // Author: Dmitry Sorokin, Google Summer of Co 27 // Author: Dmitry Sorokin, Google Summer of Code 2017 28 // Supervision: John Apostolakis, CERN 28 // Supervision: John Apostolakis, CERN 29 // ------------------------------------------- 29 // -------------------------------------------------------------------- 30 30 31 #include "G4FieldUtils.hh" 31 #include "G4FieldUtils.hh" 32 32 33 #include "G4SystemOfUnits.hh" 33 #include "G4SystemOfUnits.hh" 34 #include "G4PhysicalConstants.hh" 34 #include "G4PhysicalConstants.hh" 35 35 36 namespace field_utils { 36 namespace field_utils { 37 37 38 G4double absoluteError(const G4double y[], 38 G4double absoluteError(const G4double y[], 39 const G4double yError[] 39 const G4double yError[], 40 G4double hstep) 40 G4double hstep) 41 { 41 { 42 const G4double momentum2 = getValue2(y, Val 42 const G4double momentum2 = getValue2(y, Value3D::Momentum); 43 const G4double invMomentum2 = 1.0 / momentu 43 const G4double invMomentum2 = 1.0 / momentum2; 44 const G4double positionError2 = getValue2(y 44 const G4double positionError2 = getValue2(yError, Value3D::Position); 45 const G4double momentumError2 = getValue2(y 45 const G4double momentumError2 = getValue2(yError, Value3D::Momentum); 46 const G4double relativeMomentumError2 = mom 46 const G4double relativeMomentumError2 = momentumError2 * invMomentum2; 47 47 48 return std::max(std::sqrt(positionError2), 48 return std::max(std::sqrt(positionError2), 49 std::sqrt(relativeMomentumE 49 std::sqrt(relativeMomentumError2) * hstep); 50 } 50 } 51 51 52 G4double relativeError2(const G4double y[], 52 G4double relativeError2(const G4double y[], 53 const G4double yerr[], 53 const G4double yerr[], 54 G4double h, 54 G4double h, 55 G4double eps_rel 55 G4double eps_rel_max) 56 { 56 { 57 G4double errmax_sq; 57 G4double errmax_sq; 58 58 59 G4double inv_eps_vel_sq = 1.0 / (eps_rel_ma 59 G4double inv_eps_vel_sq = 1.0 / (eps_rel_max * eps_rel_max); 60 G4double errvel_sq = 0.0; // square of m 60 G4double errvel_sq = 0.0; // square of momentum vector difference 61 61 62 G4double eps_pos = eps_rel_max * h; 62 G4double eps_pos = eps_rel_max * h; 63 G4double inv_eps_pos_sq = 1.0 / (eps_pos * 63 G4double inv_eps_pos_sq = 1.0 / (eps_pos * eps_pos); 64 64 65 // Evaluate accuracy 65 // Evaluate accuracy 66 // 66 // 67 G4double errpos_sq = getValue2(yerr, Value3 67 G4double errpos_sq = getValue2(yerr, Value3D::Position); 68 errpos_sq *= inv_eps_pos_sq; // Scale relat 68 errpos_sq *= inv_eps_pos_sq; // Scale relative to required tolerance 69 69 70 // Accuracy for momentum 70 // Accuracy for momentum 71 // 71 // 72 G4double magvel_sq = getValue2(y, Value3D:: 72 G4double magvel_sq = getValue2(y, Value3D::Momentum); 73 G4double sumerr_sq = getValue2(yerr, Value3 73 G4double sumerr_sq = getValue2(yerr, Value3D::Momentum); 74 if (magvel_sq > 0.0) 74 if (magvel_sq > 0.0) 75 { 75 { 76 errvel_sq = sumerr_sq / magvel_sq; 76 errvel_sq = sumerr_sq / magvel_sq; 77 } 77 } 78 else 78 else 79 { 79 { 80 G4Exception("field_utils::relativeError" 80 G4Exception("field_utils::relativeError","Field001", 81 JustWarning, "found case of 81 JustWarning, "found case of zero momentum"); 82 errvel_sq = sumerr_sq; 82 errvel_sq = sumerr_sq; 83 } 83 } 84 errvel_sq *= inv_eps_vel_sq; 84 errvel_sq *= inv_eps_vel_sq; 85 errmax_sq = std::max(errpos_sq, errvel_sq); 85 errmax_sq = std::max(errpos_sq, errvel_sq); 86 86 87 return errmax_sq; 87 return errmax_sq; 88 } 88 } 89 89 90 G4double relativeError(const G4double y[], 90 G4double relativeError(const G4double y[], 91 const G4double yError[] 91 const G4double yError[], 92 const G4double h, 92 const G4double h, 93 const G4double errorTol 93 const G4double errorTolerance) 94 { 94 { 95 return std::sqrt(relativeError2(y, yError, 95 return std::sqrt(relativeError2(y, yError, h, errorTolerance)); 96 } 96 } 97 97 98 void copy(G4double dst[], const G4double src[] << 98 void copy(G4double dst[], const G4double src[], size_t size) 99 { 99 { 100 std::memcpy(dst, src, sizeof(G4double) * si 100 std::memcpy(dst, src, sizeof(G4double) * size); 101 } 101 } 102 102 103 103 104 G4double inverseCurvatureRadius(G4double parti 104 G4double inverseCurvatureRadius(G4double particleCharge, 105 G4double momen 105 G4double momentum, 106 G4double BFiel 106 G4double BField) 107 { 107 { 108 return -c_light * particleCharge * BField / 108 return -c_light * particleCharge * BField / momentum; 109 } 109 } 110 110 111 } // field_utils 111 } // field_utils 112 112 113 113