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