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