| 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 <<
26 // 25 //
27 // Author: Dmitry Sorokin, Google Summer of Co << 26 // $Id: $
28 // Supervision: John Apostolakis, CERN << 27 //
29 // ------------------------------------------- << 28 //
>> 29 // Implementation by Dmitry Sorokin - GSoC 2017
>> 30 // Work supported by Google as part of Google Summer of Code 2017.
>> 31 // Supervision / code review: John Apostolakis
30 32
31 #include "G4FieldUtils.hh" <<
32 33
33 #include "G4SystemOfUnits.hh" << 34 #include "G4FieldUtils.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 relativeError(
39 const G4double yError[] << 39 const G4double y[],
40 G4double hstep) << 40 const G4double yError[],
41 { << 41 const G4double h,
42 const G4double momentum2 = getValue2(y, Val << 42 const G4double errorTolerance)
43 const G4double invMomentum2 = 1.0 / momentu <<
44 const G4double positionError2 = getValue2(y <<
45 const G4double momentumError2 = getValue2(y <<
46 const G4double relativeMomentumError2 = mom <<
47 <<
48 return std::max(std::sqrt(positionError2), <<
49 std::sqrt(relativeMomentumE <<
50 } <<
51 <<
52 G4double relativeError2(const G4double y[], <<
53 const G4double yerr[], <<
54 G4double h, <<
55 G4double eps_rel <<
56 { <<
57 G4double errmax_sq; <<
58 <<
59 G4double inv_eps_vel_sq = 1.0 / (eps_rel_ma <<
60 G4double errvel_sq = 0.0; // square of m <<
61 <<
62 G4double eps_pos = eps_rel_max * h; <<
63 G4double inv_eps_pos_sq = 1.0 / (eps_pos * <<
64 <<
65 // Evaluate accuracy <<
66 // <<
67 G4double errpos_sq = getValue2(yerr, Value3 <<
68 errpos_sq *= inv_eps_pos_sq; // Scale relat <<
69 <<
70 // Accuracy for momentum <<
71 // <<
72 G4double magvel_sq = getValue2(y, Value3D:: <<
73 G4double sumerr_sq = getValue2(yerr, Value3 <<
74 if (magvel_sq > 0.0) <<
75 { <<
76 errvel_sq = sumerr_sq / magvel_sq; <<
77 } <<
78 else <<
79 { <<
80 G4Exception("field_utils::relativeError" <<
81 JustWarning, "found case of <<
82 errvel_sq = sumerr_sq; <<
83 } <<
84 errvel_sq *= inv_eps_vel_sq; <<
85 errmax_sq = std::max(errpos_sq, errvel_sq); <<
86 <<
87 return errmax_sq; <<
88 } <<
89 <<
90 G4double relativeError(const G4double y[], <<
91 const G4double yError[] <<
92 const G4double h, <<
93 const G4double errorTol <<
94 { 43 {
95 return std::sqrt(relativeError2(y, yError, << 44 // Accuracy for position
96 } << 45 G4double error2 = getValue2(yError, Value3D::Position) / sqr(h);
97 46
98 void copy(G4double dst[], const G4double src[] << 47 // Accuracy for momentum
99 { << 48 const G4double momentum2 = getValue2(y, Value3D::Momentum);
100 std::memcpy(dst, src, sizeof(G4double) * si << 49 if (momentum2 > 0) {
101 } << 50 const G4double momentumError2 =
102 << 51 getValue2(yError, Value3D::Momentum) / momentum2;
103 << 52 error2 = std::max(error2, momentumError2);
104 G4double inverseCurvatureRadius(G4double parti << 53 } else {
105 G4double momen << 54 G4Exception("field_utils::relativeError","Field001",
106 G4double BFiel << 55 JustWarning, "found case of zero momentum");
107 { << 56 }
108 return -c_light * particleCharge * BField / << 57 #if 0
>> 58 // Accuracy for spin
>> 59 const G4double spin2 = getValue2(y, Value3D::Spin);
>> 60 if (spin2 > 0) {
>> 61 const G4double spinError2 = getValue2(yError, Value3D::Spin) / spin2;
>> 62 error2 = std::max(error2, spinError2);
>> 63 }
>> 64 #endif
>> 65 return std::sqrt(error2) / errorTolerance;
109 } 66 }
110 67
111 } // field_utils 68 } // field_utils
112 69
113 70