Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/global/management/include/G4Pow.hh

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /global/management/include/G4Pow.hh (Version 11.3.0) and /global/management/include/G4Pow.hh (Version 9.6)


  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 // G4Pow                                       <<  26 // $Id$
                                                   >>  27 //
                                                   >>  28 //
                                                   >>  29 // -------------------------------------------------------------------
                                                   >>  30 //
                                                   >>  31 // Class G4Pow
 27 //                                                 32 //
 28 // Class description:                              33 // Class description:
 29 //                                                 34 //
 30 // Utility singleton class for the fast comput     35 // Utility singleton class for the fast computation of log and pow
 31 // functions. Integer argument should be in th <<  36 // functions. Integer argument should in the interval 0-512, no
 32 // check is performed inside these methods for     37 // check is performed inside these methods for performance reasons.
 33 // For factorial integer argument should be in     38 // For factorial integer argument should be in the interval 0-170
 34 // Computations with double arguments are fast     39 // Computations with double arguments are fast for the interval
 35 // 0.002-511.5 for all functions except expone <<  40 // 0.5-255.5, standard library is used in the opposite case
 36 // for the interval 0-84.4, standard library i <<  41 
                                                   >>  42 // Author: Vladimir Ivanchenko
                                                   >>  43 //
                                                   >>  44 // Creation date: 23.05.2009
                                                   >>  45 // -------------------------------------------------------------------
 37                                                    46 
 38 // Author: Vladimir Ivanchenko, 23.05.2009     <<  47 #ifndef G4Pow_h
 39 // ------------------------------------------- <<  48 #define G4Pow_h 1
 40 #ifndef G4Pow_hh                               << 
 41 #define G4Pow_hh 1                             << 
 42                                                    49 
 43 #include "G4DataVector.hh"                     << 
 44 #include "G4Exp.hh"                            << 
 45 #include "G4Log.hh"                            << 
 46 #include "globals.hh"                              50 #include "globals.hh"
                                                   >>  51 #include "G4DataVector.hh"
 47                                                    52 
 48 class G4Pow                                        53 class G4Pow
 49 {                                                  54 {
 50  public:                                       << 
 51   static G4Pow* GetInstance();                 << 
 52   ~G4Pow() = default;                          << 
 53                                                << 
 54   // Fast computation of Z^1/3                 << 
 55   //                                           << 
 56   inline G4double Z13(G4int Z) const;          << 
 57   G4double A13(G4double A) const;              << 
 58                                                << 
 59   // Fast computation of Z^2/3                 << 
 60   //                                           << 
 61   inline G4double Z23(G4int Z) const;          << 
 62   inline G4double A23(G4double A) const;       << 
 63                                                << 
 64   // Fast computation of log(Z)                << 
 65   //                                           << 
 66   inline G4double logZ(G4int Z) const;         << 
 67   inline G4double logA(G4double A) const;      << 
 68   inline G4double logX(G4double x) const;      << 
 69                                                << 
 70   // Fast computation of log10(Z)              << 
 71   //                                           << 
 72   inline G4double log10Z(G4int Z) const;       << 
 73   inline G4double log10A(G4double A) const;    << 
 74                                                << 
 75   // Fast computation of exp(X)                << 
 76   //                                           << 
 77   inline G4double expA(G4double A) const;      << 
 78                                                << 
 79   // Fast computation of pow(Z,X)              << 
 80   //                                           << 
 81   inline G4double powZ(G4int Z, G4double y) co << 
 82   inline G4double powA(G4double A, G4double y) << 
 83   G4double powN(G4double x, G4int n) const;    << 
 84                                                << 
 85   // Fast factorial                            << 
 86   //                                           << 
 87   inline G4double factorial(G4int Z) const;    << 
 88   inline G4double logfactorial(G4int Z) const; << 
 89                                                << 
 90  private:                                      << 
 91   G4Pow();                                     << 
 92                                                << 
 93   G4double A13Low(const G4double, const G4bool << 
 94   G4double A13High(const G4double, const G4boo << 
 95                                                << 
 96   inline G4double logBase(G4double x) const;   << 
 97                                                << 
 98   static G4Pow* fpInstance;                    << 
 99                                                << 
100   const G4double onethird = 1.0 / 3.0;         << 
101   const G4int max2        = 5;                 << 
102                                                << 
103   G4double maxA;                               << 
104   G4double maxLowA;                            << 
105   G4double maxA2;                              << 
106   G4double maxAexp;                            << 
107                                                << 
108   G4DataVector ener;                           << 
109   G4DataVector logen;                          << 
110   G4DataVector pz13;                           << 
111   G4DataVector lowa13;                         << 
112   G4DataVector lz;                             << 
113   G4DataVector lz2;                            << 
114   G4DataVector fexp;                           << 
115   G4DataVector fact;                           << 
116   G4DataVector logfact;                        << 
117 };                                             << 
118                                                    55 
119 // -----------------------------               <<  56   public:
120 // Inline methods implementation               << 
121 // -----------------------------               << 
122                                                    57 
123 inline G4double G4Pow::Z13(G4int Z) const { re <<  58     static G4Pow* GetInstance();
124                                                    59 
125 inline G4double G4Pow::Z23(G4int Z) const      <<  60     // Fast computation of Z^1/3
126 {                                              <<  61     //
127   G4double x = Z13(Z);                         <<  62     inline G4double Z13(G4int Z);
128   return x * x;                                <<  63     inline G4double A13(G4double A);
129 }                                              <<  64 
                                                   >>  65     // Fast computation of Z^2/3
                                                   >>  66     //
                                                   >>  67     inline G4double Z23(G4int Z);
                                                   >>  68     inline G4double A23(G4double A);
                                                   >>  69 
                                                   >>  70     // Fast computation of log(Z)
                                                   >>  71     //
                                                   >>  72     inline G4double logZ(G4int Z);
                                                   >>  73     inline G4double logA(G4double A);
                                                   >>  74 
                                                   >>  75     // Fast computation of log10(Z)
                                                   >>  76     //
                                                   >>  77     inline G4double log10Z(G4int Z);
                                                   >>  78     inline G4double log10A(G4double A);
                                                   >>  79 
                                                   >>  80     // Fast computation of pow(Z,X)
                                                   >>  81     //
                                                   >>  82     inline G4double powZ(G4int Z, G4double y);
                                                   >>  83     inline G4double powA(G4double A, G4double y);
                                                   >>  84            G4double powN(G4double x, G4int n);
                                                   >>  85 
                                                   >>  86     // Fast factorial
                                                   >>  87     //
                                                   >>  88     inline G4double factorial(G4int Z);
                                                   >>  89     inline G4double logfactorial(G4int Z);
130                                                    90 
131 inline G4double G4Pow::A23(G4double A) const   <<  91   private:
                                                   >>  92 
                                                   >>  93     G4Pow();
                                                   >>  94    ~G4Pow();
                                                   >>  95 
                                                   >>  96   private:
                                                   >>  97 
                                                   >>  98     static G4Pow* fInstance;
                                                   >>  99 
                                                   >> 100     G4double onethird;
                                                   >> 101 
                                                   >> 102     G4DataVector pz13;
                                                   >> 103     G4DataVector lz;
                                                   >> 104     G4DataVector fact;
                                                   >> 105     G4DataVector logfact;
                                                   >> 106 };
                                                   >> 107 
                                                   >> 108 // -------------------------------------------------------------------
                                                   >> 109 
                                                   >> 110 inline G4double G4Pow::Z13(G4int Z)
132 {                                                 111 {
133   G4double x = A13(A);                         << 112   return pz13[Z];
134   return x * x;                                << 
135 }                                                 113 }
136                                                   114 
137 inline G4double G4Pow::logZ(G4int Z) const { r << 115 inline G4double G4Pow::A13(G4double A)
138                                                << 
139 inline G4double G4Pow::logBase(G4double a) con << 
140 {                                                 116 {
                                                   >> 117   const G4double minA = 0.5000001; 
                                                   >> 118   const G4double maxA = 255.5; 
                                                   >> 119 
141   G4double res;                                   120   G4double res;
142   if(a <= maxA2)                               << 121   if((A >= minA) && (A <= maxA))
143   {                                               122   {
144     G4int i = G4int(max2 * (a - 1) + 0.5);     << 123     G4int i = G4int(A + 0.5);
145     if(i > max2)                               << 124     G4double x = (1.0 - A/G4double(i))*onethird;
146     {                                          << 125     res = pz13[i]*(1.0 + x - x*x*(1.0 - 1.66666666*x));
147       i = max2;                                << 
148     }                                          << 
149     G4double x = a / (G4double(i) / max2 + 1)  << 
150     res        = lz2[i] + x * (1.0 - (0.5 - on << 
151   }                                            << 
152   else if(a <= maxA)                           << 
153   {                                            << 
154     G4int i    = G4int(a + 0.5);               << 
155     G4double x = a / G4double(i) - 1;          << 
156     res        = lz[i] + x * (1.0 - (0.5 - one << 
157   }                                               126   }
158   else                                            127   else
159   {                                               128   {
160     res = G4Log(a);                            << 129     res = std::pow(A, onethird); 
161   }                                               130   }
162   return res;                                     131   return res;
163 }                                                 132 }
164                                                   133 
165 inline G4double G4Pow::logA(G4double A) const  << 134 inline G4double G4Pow::Z23(G4int Z)
                                                   >> 135 {
                                                   >> 136   G4double x = Z13(Z);
                                                   >> 137   return x*x;
                                                   >> 138 }
                                                   >> 139 
                                                   >> 140 inline G4double G4Pow::A23(G4double A)
166 {                                                 141 {
167   return (1.0 <= A ? logBase(A) : -logBase(1.  << 142   G4double x = A13(A);
                                                   >> 143   return x*x;
168 }                                                 144 }
169                                                   145 
170 inline G4double G4Pow::logX(G4double x) const  << 146 inline G4double G4Pow::logZ(G4int Z)
171 {                                                 147 {
172   G4double res = 0.0;                          << 148   return lz[Z];
173   G4double a   = (1.0 <= x) ? x : 1.0 / x;     << 149 }
174                                                   150 
175   if(a <= maxA)                                << 151 inline G4double G4Pow::logA(G4double A)
176   {                                            << 152 {
177     res = logBase(a);                          << 153   const G4double minA = 0.5000001; 
178   }                                            << 154   const G4double maxA = 255.5; 
179   else if(a <= ener[2])                        << 155 
180   {                                            << 156   G4double res;
181     res = logen[1] + logBase(a / ener[1]);     << 157   if((A >= minA) && (A <= maxA))
182   }                                            << 
183   else if(a <= ener[3])                        << 
184   {                                               158   {
185     res = logen[2] + logBase(a / ener[2]);     << 159     G4int i = G4int(A + 0.5);
                                                   >> 160     G4double x = 1.0 - A/G4double(i);
                                                   >> 161     res = lz[i] + x*(1.0 - 0.5*x + onethird*x*x);
186   }                                               162   }
187   else                                            163   else
188   {                                               164   {
189     res = G4Log(a);                            << 165     res = std::log(A);
190   }                                            << 
191                                                << 
192   if(1.0 > x)                                  << 
193   {                                            << 
194     res = -res;                                << 
195   }                                               166   }
196   return res;                                     167   return res;
197 }                                                 168 }
198                                                   169 
199 inline G4double G4Pow::log10Z(G4int Z) const { << 170 inline G4double G4Pow::log10Z(G4int Z)
200                                                << 171 {
201 inline G4double G4Pow::log10A(G4double A) cons << 172   return lz[Z]/lz[10];
                                                   >> 173 }
202                                                   174 
203 inline G4double G4Pow::expA(G4double A) const  << 175 inline G4double G4Pow::log10A(G4double A)
204 {                                                 176 {
205   G4double res;                                << 177   return logA(A)/lz[10];
206   G4double a = (0.0 <= A) ? A : -A;            << 178 }
207                                                   179 
208   if(a <= maxAexp)                             << 180 inline G4double G4Pow::powZ(G4int Z, G4double y)
209   {                                            << 181 {
210     G4int i    = G4int(2 * a + 0.5);           << 182   return std::exp(y*lz[Z]);
211     G4double x = a - i * 0.5;                  << 
212     res        = fexp[i] * (1.0 + x * (1.0 + 0 << 
213   }                                            << 
214   else                                         << 
215   {                                            << 
216     res = G4Exp(a);                            << 
217   }                                            << 
218   if(0.0 > A)                                  << 
219   {                                            << 
220     res = 1.0 / res;                           << 
221   }                                            << 
222   return res;                                  << 
223 }                                                 183 }
224                                                   184 
225 inline G4double G4Pow::powZ(G4int Z, G4double  << 185 inline G4double G4Pow::powA(G4double A, G4double y)
226 {                                                 186 {
227   return expA(y * lz[Z]);                      << 187   return std::exp(y*logA(A));
228 }                                                 188 }
229                                                   189 
230 inline G4double G4Pow::powA(G4double A, G4doub << 190 inline G4double G4Pow::factorial(G4int Z)
231 {                                                 191 {
232   return (0.0 == A ? 0.0 : expA(y * logX(A))); << 192   return fact[Z];
233 }                                                 193 }
234                                                   194 
235 inline G4double G4Pow::factorial(G4int Z) cons << 195 inline G4double G4Pow::logfactorial(G4int Z)
                                                   >> 196 {
                                                   >> 197   return logfact[Z];
                                                   >> 198 }
236                                                   199 
237 inline G4double G4Pow::logfactorial(G4int Z) c << 200 // -------------------------------------------------------------------
238                                                   201 
239 #endif                                            202 #endif
                                                   >> 203 
240                                                   204