| 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 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France 27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium 28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France 29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H 30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France 31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France
33 // 33 //
34 #define INCLXX_IN_GEANT4_MODE 1 34 #define INCLXX_IN_GEANT4_MODE 1
35 35
36 #include "globals.hh" 36 #include "globals.hh"
37 37
38 #include "G4INCLGlobals.hh" 38 #include "G4INCLGlobals.hh"
39 #include "G4INCLParticle.hh" 39 #include "G4INCLParticle.hh"
40 #ifdef HAVE_WORDEXP 40 #ifdef HAVE_WORDEXP
41 #include <wordexp.h> 41 #include <wordexp.h>
42 #endif 42 #endif
43 43
44 namespace G4INCL { 44 namespace G4INCL {
45 namespace Math { 45 namespace Math {
46 46
47 namespace { 47 namespace {
48 48
49 // constants for the Gaussian CDF approx 49 // constants for the Gaussian CDF approximation
50 const G4double gcdfa1 = 0.254829592; 50 const G4double gcdfa1 = 0.254829592;
51 const G4double gcdfa2 = -0.284496736; 51 const G4double gcdfa2 = -0.284496736;
52 const G4double gcdfa3 = 1.421413741; 52 const G4double gcdfa3 = 1.421413741;
53 const G4double gcdfa4 = -1.453152027; 53 const G4double gcdfa4 = -1.453152027;
54 const G4double gcdfa5 = 1.061405429; 54 const G4double gcdfa5 = 1.061405429;
55 const G4double gcdfp = 0.3275911; 55 const G4double gcdfp = 0.3275911;
56 56
57 // constants for the inverse Gaussian CD 57 // constants for the inverse Gaussian CDF approximation
58 const G4double igcdfc1 = 2.515517; 58 const G4double igcdfc1 = 2.515517;
59 const G4double igcdfc2 = 0.802853; 59 const G4double igcdfc2 = 0.802853;
60 const G4double igcdfc3 = 0.010328; 60 const G4double igcdfc3 = 0.010328;
61 const G4double igcdfd1 = 1.432788; 61 const G4double igcdfd1 = 1.432788;
62 const G4double igcdfd2 = 0.189269; 62 const G4double igcdfd2 = 0.189269;
63 const G4double igcdfd3 = 0.001308; 63 const G4double igcdfd3 = 0.001308;
64 64
65 G4double inverseGaussianCDFRational(cons 65 G4double inverseGaussianCDFRational(const G4double t) {
66 // Abramowitz and Stegun formula 26.2. 66 // Abramowitz and Stegun formula 26.2.23.
67 // The absolute value of the error sho 67 // The absolute value of the error should be less than 4.5 e-4.
68 return t - ((igcdfc3*t + igcdfc2)*t + 68 return t - ((igcdfc3*t + igcdfc2)*t + igcdfc1) /
69 (((igcdfd3*t + igcdfd2)*t + igcdfd1) 69 (((igcdfd3*t + igcdfd2)*t + igcdfd1)*t + 1.0);
70 } 70 }
71 71
72 } 72 }
73 73
74 G4double gaussianCDF(const G4double x) 74 G4double gaussianCDF(const G4double x)
75 { 75 {
76 // Save the sign of x 76 // Save the sign of x
77 const G4double sgn = sign(x); 77 const G4double sgn = sign(x);
78 const G4double z = std::fabs(x) * oneOve 78 const G4double z = std::fabs(x) * oneOverSqrtTwo;
79 79
80 // A&S formula 7.1.26 80 // A&S formula 7.1.26
81 G4double t = 1.0/(1.0 + gcdfp*z); 81 G4double t = 1.0/(1.0 + gcdfp*z);
82 G4double y = 1.0 - (((((gcdfa5*t + gcdfa 82 G4double y = 1.0 - (((((gcdfa5*t + gcdfa4)*t) + gcdfa3)*t + gcdfa2)*t + gcdfa1)*t*std::exp(-z*z);
83 83
84 return 0.5*(1.0 + sgn*y); 84 return 0.5*(1.0 + sgn*y);
85 } 85 }
86 86
87 G4double gaussianCDF(const G4double x, con 87 G4double gaussianCDF(const G4double x, const G4double x0, const G4double sigma) {
88 return gaussianCDF((x-x0)/sigma); 88 return gaussianCDF((x-x0)/sigma);
89 } 89 }
90 90
91 G4double inverseGaussianCDF(G4double x) { 91 G4double inverseGaussianCDF(G4double x) {
92 if (x < 0.5) 92 if (x < 0.5)
93 return -inverseGaussianCDFRational( st 93 return -inverseGaussianCDFRational( std::sqrt(-2.0*std::log(x)) );
94 else 94 else
95 return inverseGaussianCDFRational( std 95 return inverseGaussianCDFRational( std::sqrt(-2.0*std::log(1.-x)) );
96 } 96 }
97 97
98 G4double arcSin(const G4double x) { 98 G4double arcSin(const G4double x) {
99 // assert(x>-1.000001 && x<1.000001); 99 // assert(x>-1.000001 && x<1.000001);
100 return ((x > 1.) ? 0. : ((x<-1.) ? pi : 100 return ((x > 1.) ? 0. : ((x<-1.) ? pi : std::asin(x)));
101 } 101 }
102 102
103 G4double arcCos(const G4double x) { 103 G4double arcCos(const G4double x) {
104 // assert(x>-1.000001 && x<1.000001); 104 // assert(x>-1.000001 && x<1.000001);
105 return ((x > 1.) ? 0. : ((x<-1.) ? pi : 105 return ((x > 1.) ? 0. : ((x<-1.) ? pi : std::acos(x)));
106 } 106 }
107 } 107 }
108 108
109 namespace ParticleConfig { 109 namespace ParticleConfig {
110 G4bool isPair(Particle const * const p1, P 110 G4bool isPair(Particle const * const p1, Particle const * const p2, ParticleType t1, ParticleType t2) {
111 return ((p1->getType() == t1 && p2->getT 111 return ((p1->getType() == t1 && p2->getType() == t2) ||
112 (p1->getType() == t2 && p2->getType() 112 (p1->getType() == t2 && p2->getType() == t1));
113 } 113 }
114 } 114 }
115 115
116 #ifndef INCLXX_IN_GEANT4_MODE 116 #ifndef INCLXX_IN_GEANT4_MODE
117 namespace String { 117 namespace String {
118 void wrap(std::string &str, const size_t l 118 void wrap(std::string &str, const size_t lineLength, const std::string &separators) {
119 const size_t len = str.size(); 119 const size_t len = str.size();
120 size_t startPos = 0; 120 size_t startPos = 0;
121 while(len-startPos > lineLength) { /* Lo 121 while(len-startPos > lineLength) { /* Loop checking, 10.07.2015, D.Mancusi */
122 const size_t nextNewline = str.find('\ 122 const size_t nextNewline = str.find('\n', startPos);
123 if(nextNewline!=std::string::npos && n 123 if(nextNewline!=std::string::npos && nextNewline-startPos<=lineLength)
124 startPos = nextNewline+1; 124 startPos = nextNewline+1;
125 else { 125 else {
126 size_t lastSeparator = str.find_last 126 size_t lastSeparator = str.find_last_of(separators, startPos+lineLength);
127 if(lastSeparator!=std::string::npos) 127 if(lastSeparator!=std::string::npos)
128 str[lastSeparator] = '\n'; 128 str[lastSeparator] = '\n';
129 startPos = lastSeparator+1; 129 startPos = lastSeparator+1;
130 } 130 }
131 } 131 }
132 } 132 }
133 133
134 void replaceAll(std::string &str, const st 134 void replaceAll(std::string &str, const std::string &from, const std::string &to, const size_t maxPosition) {
135 if(from.empty()) 135 if(from.empty())
136 return; 136 return;
137 size_t start_pos = 0; 137 size_t start_pos = 0;
138 size_t cur_max_pos = maxPosition; 138 size_t cur_max_pos = maxPosition;
139 const size_t from_len = from.length(); 139 const size_t from_len = from.length();
140 const size_t to_len = to.length(); 140 const size_t to_len = to.length();
141 while((start_pos = str.find(from, start_ 141 while((start_pos = str.find(from, start_pos)) != std::string::npos /* Loop checking, 10.07.2015, D.Mancusi */
142 && (cur_max_pos==std::string::npos 142 && (cur_max_pos==std::string::npos || start_pos<cur_max_pos)) {
143 str.replace(start_pos, from_len, to); 143 str.replace(start_pos, from_len, to);
144 start_pos += to_len; // In case 'to' c 144 start_pos += to_len; // In case 'to' contains 'from', like replacing 'x' with 'yx'
145 if(cur_max_pos!=std::string::npos) 145 if(cur_max_pos!=std::string::npos)
146 cur_max_pos += to_len - from_len; 146 cur_max_pos += to_len - from_len;
147 } 147 }
148 } 148 }
149 149
150 std::vector<std::string> tokenize(std::str 150 std::vector<std::string> tokenize(std::string const &str, const std::string &delimiters) {
151 size_t startPos = 0, endPos; 151 size_t startPos = 0, endPos;
152 std::vector<std::string> tokens; 152 std::vector<std::string> tokens;
153 do { 153 do {
154 endPos = str.find_first_of(delimiters, 154 endPos = str.find_first_of(delimiters, startPos);
155 std::string token = str.substr(startPo 155 std::string token = str.substr(startPos, endPos-startPos);
156 tokens.push_back(token); 156 tokens.push_back(token);
157 startPos = str.find_first_not_of(delim 157 startPos = str.find_first_not_of(delimiters, endPos);
158 } while(endPos!=std::string::npos); /* L 158 } while(endPos!=std::string::npos); /* Loop checking, 10.07.2015, D.Mancusi */
159 159
160 return tokens; 160 return tokens;
161 } 161 }
162 162
163 G4bool isInteger(std::string const &str) { 163 G4bool isInteger(std::string const &str) {
164 const size_t pos = str.find_first_not_of 164 const size_t pos = str.find_first_not_of("0123456789");
165 return (pos==std::string::npos); 165 return (pos==std::string::npos);
166 } 166 }
167 167
168 std::string expandPath(std::string const & 168 std::string expandPath(std::string const &path) {
169 #ifdef HAVE_WORDEXP 169 #ifdef HAVE_WORDEXP
170 wordexp_t expansionResult; 170 wordexp_t expansionResult;
171 std::string result; 171 std::string result;
172 G4int err = wordexp(path.c_str(), &expan 172 G4int err = wordexp(path.c_str(), &expansionResult, WRDE_NOCMD);
173 if(err) 173 if(err)
174 result = path; 174 result = path;
175 else 175 else
176 result = expansionResult.we_wordv[0]; 176 result = expansionResult.we_wordv[0];
177 wordfree(&expansionResult); 177 wordfree(&expansionResult);
178 return result; 178 return result;
179 #else 179 #else
180 // no-op if wordexp.h is not found 180 // no-op if wordexp.h is not found
181 return path; 181 return path;
182 #endif 182 #endif
183 } 183 }
184 } 184 }
185 185
186 #endif // INCLXX_IN_GEANT4_MODE 186 #endif // INCLXX_IN_GEANT4_MODE
187 187
188 } 188 }
189 189