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 /** \file G4INCLRootFinder.hh 37 /** \file G4INCLRootFinder.hh 39 * \brief Static root-finder algorithm. 38 * \brief Static root-finder algorithm. 40 * 39 * 41 * Provides a stateless root-finder algorithm. 40 * Provides a stateless root-finder algorithm. 42 * 41 * 43 * \date 2nd March 2011 42 * \date 2nd March 2011 44 * \author Davide Mancusi 43 * \author Davide Mancusi 45 */ 44 */ 46 45 47 #include "G4INCLRootFinder.hh" 46 #include "G4INCLRootFinder.hh" 48 #include "G4INCLGlobals.hh" 47 #include "G4INCLGlobals.hh" 49 #include "G4INCLLogger.hh" 48 #include "G4INCLLogger.hh" 50 #include <utility> 49 #include <utility> 51 #include <cmath> 50 #include <cmath> 52 51 53 namespace G4INCL { 52 namespace G4INCL { 54 53 55 namespace RootFinder { << 54 std::pair<G4double,G4double> RootFinder::solution; 56 55 57 namespace { << 56 const G4double RootFinder::toleranceY = 1.e-4; 58 57 59 /// \brief Tolerance on the y value << 58 G4bool RootFinder::solve(RootFunctor const * const f, const G4double x0) { 60 const G4double toleranceY = 1.e-4; << 59 // If we already have the solution, do nothing >> 60 const G4double y0 = (*f)(x0); >> 61 if( std::abs(y0) < toleranceY ) { >> 62 solution = std::make_pair(x0,y0); >> 63 return true; >> 64 } 61 65 62 /// \brief Maximum number of iterations << 66 // Bracket the root and set the initial values 63 const G4int maxIterations=50; << 67 std::pair<G4double,G4double> bracket = bracketRoot(f,x0); 64 << 68 G4double x1 = bracket.first; 65 /** \brief Bracket the root of the funct << 69 G4double x2 = bracket.second; 66 * << 70 // If x1>x2, it means that we could not bracket the root. Return false. 67 * Tries to find a bracketing value for << 71 if(x1>x2) { 68 * << 72 // Maybe zero is a good solution? 69 * \param f pointer to a RootFunctor << 73 G4double y_at_zero = (*f)(0.); 70 * \param x0 starting value << 74 if(std::abs(y_at_zero)<=toleranceY) { 71 * \return if the root could be brackete << 75 f->cleanUp(true); 72 * bracketing the root, as a pair. If << 76 solution = std::make_pair(0.,y_at_zero); 73 * pair with first > second. << 77 return true; 74 */ << 78 } else { 75 std::pair<G4double,G4double> bracketRoot << 79 WARN("Root-finding algorithm could not bracket the root." << std::endl); 76 G4double y0 = (*f)(x0); << 80 f->cleanUp(false); 77 << 81 return false; 78 const G4double scaleFactor = 1.5; << 79 << 80 G4double x1; << 81 if(x0!=0.) << 82 x1=scaleFactor*x0; << 83 else << 84 x1=1.; << 85 G4double y1 = (*f)(x1); << 86 << 87 if(Math::sign(y0)!=Math::sign(y1)) << 88 return std::make_pair(x0,x1); << 89 << 90 const G4double scaleFactorMinus1 = 1./ << 91 G4double oldx0, oldx1, oldy1; << 92 G4int iterations=0; << 93 do { << 94 if(iterations > maxIterations) { << 95 INCL_DEBUG("Could not bracket the << 96 return std::make_pair((G4double) 1 << 97 } << 98 << 99 oldx0=x0; << 100 oldx1=x1; << 101 oldy1=y1; << 102 << 103 x0 *= scaleFactorMinus1; << 104 x1 *= scaleFactor; << 105 y0 = (*f)(x0); << 106 y1 = (*f)(x1); << 107 iterations++; << 108 } while(Math::sign(y0)==Math::sign(y1) << 109 << 110 if(Math::sign(y1)==Math::sign(oldy1)) << 111 return std::make_pair(x0,oldx0); << 112 else << 113 return std::make_pair(oldx1,x1); << 114 } 82 } 115 << 116 } 83 } 117 84 118 Solution solve(RootFunctor const * const f << 85 G4double y1 = (*f)(x1); 119 // If we already have the solution, do n << 86 G4double y2 = (*f)(x2); 120 const G4double y0 = (*f)(x0); << 87 G4double x = x1; 121 if( std::abs(y0) < toleranceY ) { << 88 G4double y = y1; 122 return Solution(x0,y0); << 89 >> 90 /* ******************************** >> 91 * Start of the false-position loop >> 92 * ********************************/ >> 93 >> 94 // Keep track of the last updated interval end (-1=left, 1=right) >> 95 G4int lastUpdated = 0; >> 96 >> 97 for(G4int iterations=0; std::abs(y) > toleranceY; iterations++) { >> 98 >> 99 if(iterations > maxIterations) { >> 100 WARN("Root-finding algorithm did not converge." << std::endl); >> 101 f->cleanUp(false); >> 102 return false; 123 } 103 } 124 104 125 // Bracket the root and set the initial << 105 // Estimate the root position by linear interpolation 126 std::pair<G4double,G4double> bracket = b << 106 x = (y1*x2-y2*x1)/(y1-y2); 127 G4double x1 = bracket.first; << 107 128 G4double x2 = bracket.second; << 108 // Update the value of the function 129 // If x1>x2, it means that we could not << 109 y = (*f)(x); 130 if(x1>x2) { << 131 // Maybe zero is a good solution? << 132 G4double y_at_zero = (*f)(0.); << 133 if(std::abs(y_at_zero)<=toleranceY) { << 134 f->cleanUp(true); << 135 return Solution(0.,y_at_zero); << 136 } else { << 137 INCL_DEBUG("Root-finding algorithm c << 138 f->cleanUp(false); << 139 return Solution(); << 140 } << 141 } << 142 110 143 G4double y1 = (*f)(x1); << 111 // Update the bracketing interval 144 G4double y2 = (*f)(x2); << 112 if(Math::sign(y) == Math::sign(y1)) { 145 G4double x = x1; << 113 x1=x; 146 G4double y = y1; << 114 y1=y; 147 << 115 if(lastUpdated==-1) y2 *= 0.5; 148 /* ******************************** << 116 lastUpdated = -1; 149 * Start of the false-position loop << 117 } else { 150 * ********************************/ << 118 x2=x; 151 << 119 y2=y; 152 // Keep track of the last updated interv << 120 if(lastUpdated==1) y1 *= 0.5; 153 G4int lastUpdated = 0; << 121 lastUpdated = 1; 154 << 155 for(G4int iterations=0; std::abs(y) > to << 156 << 157 if(iterations > maxIterations) { << 158 INCL_DEBUG("Root-finding algorithm d << 159 f->cleanUp(false); << 160 return Solution(); << 161 } << 162 << 163 // Estimate the root position by linea << 164 x = (y1*x2-y2*x1)/(y1-y2); << 165 << 166 // Update the value of the function << 167 y = (*f)(x); << 168 << 169 // Update the bracketing interval << 170 if(Math::sign(y) == Math::sign(y1)) { << 171 x1=x; << 172 y1=y; << 173 if(lastUpdated==-1) y2 *= 0.5; << 174 lastUpdated = -1; << 175 } else { << 176 x2=x; << 177 y2=y; << 178 if(lastUpdated==1) y1 *= 0.5; << 179 lastUpdated = 1; << 180 } << 181 } 122 } >> 123 } 182 124 183 /* ****************************** << 125 /* ****************************** 184 * End of the false-position loop << 126 * End of the false-position loop 185 * ******************************/ << 127 * ******************************/ >> 128 >> 129 solution = std::make_pair(x,y); >> 130 f->cleanUp(true); >> 131 return true; >> 132 } >> 133 >> 134 std::pair<G4double,G4double> RootFinder::bracketRoot(RootFunctor const * const f, G4double x0) { >> 135 G4double y0 = (*f)(x0); >> 136 >> 137 const G4double scaleFactor = 1.5; >> 138 >> 139 G4double x1; >> 140 if(x0!=0.) >> 141 x1=scaleFactor*x0; >> 142 else >> 143 x1=1.; >> 144 G4double y1 = (*f)(x1); >> 145 >> 146 if(Math::sign(y0)!=Math::sign(y1)) >> 147 return std::make_pair(x0,x1); >> 148 >> 149 const G4double scaleFactorMinus1 = 1./scaleFactor; >> 150 G4double oldx0, oldx1, oldy1; >> 151 G4int iterations=0; >> 152 do { >> 153 if(iterations > maxIterations) { >> 154 DEBUG("Could not bracket the root." << std::endl); >> 155 return std::make_pair((G4double) 1.,(G4double) -1.); >> 156 } 186 157 187 f->cleanUp(true); << 158 oldx0=x0; 188 return Solution(x,y); << 159 oldx1=x1; 189 } << 160 oldy1=y1; >> 161 >> 162 x0 *= scaleFactorMinus1; >> 163 x1 *= scaleFactor; >> 164 y0 = (*f)(x0); >> 165 y1 = (*f)(x1); >> 166 iterations++; >> 167 } while(Math::sign(y0)==Math::sign(y1)); >> 168 >> 169 if(Math::sign(y1)==Math::sign(oldy1)) >> 170 return std::make_pair(x0,oldx0); >> 171 else >> 172 return std::make_pair(oldx1,x1); >> 173 } 190 174 191 } // namespace RootFinder << 192 } 175 } 193 176