Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/global/management/src/G4Pow.cc

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/src/G4Pow.cc (Version 11.3.0) and /global/management/src/G4Pow.cc (Version 10.4.p2)


  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 class implemenation                   <<  26 // $Id: G4Pow.cc 93311 2015-10-16 10:16:37Z gcosmo $
 27 //                                                 27 //
 28 // Author: Vladimir Ivanchenko, 23.05.2009     <<  28 // -------------------------------------------------------------------
 29 // Revisions:                                  <<  29 //
 30 // - 23.03.2018, M.Novak: increased accuracy o <<  30 // GEANT4 Class file
 31 //                        [1/4,4] interval by  <<  31 //
 32 // ------------------------------------------- <<  32 //
                                                   >>  33 // File name:     G4Pow
                                                   >>  34 //
                                                   >>  35 // Author:        Vladimir Ivanchenko
                                                   >>  36 //
                                                   >>  37 // Creation date: 23.05.2009
                                                   >>  38 //
                                                   >>  39 // Modifications:
                                                   >>  40 // 08.01.2011 V.Ivanchenko extended maxZ from 256 to 512
                                                   >>  41 // 02.05.2013 V.Ivanchenko added expA and logX methods, 
                                                   >>  42 //            revised A13, logA, powZ, powA to improved accuracy
                                                   >>  43 //
                                                   >>  44 // -------------------------------------------------------------------
 33                                                    45 
 34 #include "G4Pow.hh"                                46 #include "G4Pow.hh"
                                                   >>  47 #ifdef G4MULTITHREADED
 35 #include "G4Threading.hh"                          48 #include "G4Threading.hh"
                                                   >>  49 #endif
 36                                                    50 
 37 G4Pow* G4Pow::fpInstance = nullptr;            <<  51 G4Pow* G4Pow::fpInstance = 0;
 38                                                    52 
 39 // -------------------------------------------     53 // -------------------------------------------------------------------
 40                                                    54 
 41 G4Pow* G4Pow::GetInstance()                        55 G4Pow* G4Pow::GetInstance()
 42 {                                                  56 {
 43   if(fpInstance == nullptr)                    <<  57   if (fpInstance == 0)
 44   {                                                58   {
 45     static G4Pow geant4pow;                        59     static G4Pow geant4pow;
 46     fpInstance = &geant4pow;                       60     fpInstance = &geant4pow;
 47   }                                                61   }
 48   return fpInstance;                               62   return fpInstance;
 49 }                                                  63 }
 50                                                    64 
 51 // -------------------------------------------     65 // -------------------------------------------------------------------
 52                                                    66 
 53 G4Pow::G4Pow()                                     67 G4Pow::G4Pow()
 54 {                                              <<  68   : onethird(1.0/3.0), max2(5)
                                                   >>  69 {  
 55 #ifdef G4MULTITHREADED                             70 #ifdef G4MULTITHREADED
 56   if(G4Threading::IsWorkerThread())                71   if(G4Threading::IsWorkerThread())
 57   {                                            <<  72   { 
 58     G4Exception("G4Pow::G4Pow()", "InvalidSetu <<  73     G4Exception ("G4Pow::G4Pow()", "InvalidSetup", FatalException, 
 59                 "Attempt to instantiate G4Pow  <<  74                  "Attempt to instantiate G4Pow in worker thread!");
 60   }                                                75   }
 61 #endif                                             76 #endif
 62   const G4int maxZ     = 512;                  <<  77   const G4int maxZ = 512; 
 63   const G4int maxZfact = 170;                  <<  78   const G4int maxZfact = 170; 
 64   const G4int numLowA  = 17;                   << 
 65                                                    79 
 66   maxA    = -0.6 + maxZ;                           80   maxA    = -0.6 + maxZ;
 67   maxLowA = 4.0;                               <<  81   maxA2   = 1.25 + max2*0.2;
 68   maxA2   = 1.25 + max2 * 0.2;                 <<  82   maxAexp = -0.76+ maxZfact*0.5;
 69   maxAexp = -0.76 + maxZfact * 0.5;            <<  83 
 70                                                <<  84   ener.resize(max2+1,1.0);
 71   ener.resize(max2 + 1, 1.0);                  <<  85   logen.resize(max2+1,0.0);
 72   logen.resize(max2 + 1, 0.0);                 <<  86   lz2.resize(max2+1,0.0);
 73   lz2.resize(max2 + 1, 0.0);                   <<  87   pz13.resize(maxZ,0.0);
 74   pz13.resize(maxZ, 0.0);                      <<  88   lz.resize(maxZ,0.0);
 75   lowa13.resize(numLowA, 0.0);                 <<  89   fexp.resize(maxZfact,0.0);
 76   lz.resize(maxZ, 0.0);                        <<  90   fact.resize(maxZfact,0.0);
 77   fexp.resize(maxZfact, 0.0);                  <<  91   logfact.resize(maxZ,0.0);
 78   fact.resize(maxZfact, 0.0);                  << 
 79   logfact.resize(maxZ, 0.0);                   << 
 80                                                    92 
 81   G4double f    = 1.0;                         <<  93   G4double f = 1.0;
 82   G4double logf = 0.0;                             94   G4double logf = 0.0;
 83   fact[0]       = 1.0;                         <<  95   fact[0] = 1.0;
 84   fexp[0]       = 1.0;                         <<  96   fexp[0] = 1.0;
 85                                                    97 
 86   for(G4int i = 1; i <= max2; ++i)             <<  98   for(G4int i=1; i<=max2; ++i)
 87   {                                                99   {
 88     ener[i]  = powN(500., i);                  << 100     ener[i] = powN(500.,i); 
 89     logen[i] = G4Log(ener[i]);                 << 101     logen[i]= G4Log(ener[i]); 
 90     lz2[i]   = G4Log(1.0 + i * 0.2);           << 102     lz2[i]  = G4Log(1.0 + i*0.2);
 91   }                                               103   }
 92                                                   104 
 93   for(G4int i = 1; i < maxZ; ++i)              << 105   for(G4int i=1; i<maxZ; ++i)
 94   {                                               106   {
 95     auto  x = G4double(i);                     << 107     G4double x  = G4double(i);
 96     pz13[i]    = std::pow(x, onethird);        << 108     pz13[i] = std::pow(x,onethird);
 97     lz[i]      = G4Log(x);                     << 109     lz[i]   = G4Log(x);
 98     if(i < maxZfact)                              110     if(i < maxZfact)
 99     {                                          << 111     { 
100       f *= x;                                  << 112       f *= x; 
101       fact[i] = f;                                113       fact[i] = f;
102       fexp[i] = G4Exp(0.5 * i);                << 114       fexp[i] = G4Exp(0.5*i);
103     }                                             115     }
104     logf += lz[i];                                116     logf += lz[i];
105     logfact[i] = logf;                            117     logfact[i] = logf;
106   }                                               118   }
107                                                << 
108   for(G4int i = 4; i < numLowA; ++i)           << 
109   {                                            << 
110     lowa13[i] = std::pow(0.25 * i, onethird);  << 
111   }                                            << 
112 }                                                 119 }
113                                                   120 
114 // -------------------------------------------    121 // -------------------------------------------------------------------
115                                                   122 
116 G4double G4Pow::A13(G4double A) const          << 123 G4Pow::~G4Pow()
117 {                                              << 124 {}
118   G4double res = 0.;                           << 
119   if(A > 0.)                                   << 
120   {                                            << 
121     const bool invert = (A < 1.);              << 
122     const G4double a  = invert ? 1. / A : A;   << 
123     res               = (a < maxLowA) ? A13Low << 
124   }                                            << 
125   return res;                                  << 
126 }                                              << 
127                                                << 
128 // ------------------------------------------- << 
129                                                << 
130 G4double G4Pow::A13High(const G4double a, cons << 
131 {                                              << 
132   G4double res;                                << 
133   if(a < maxA)                                 << 
134   {                                            << 
135     const auto  i    = static_cast<G4int>(a +  << 
136     const G4double x = (a / i - 1.) * onethird << 
137     res              = pz13[i] * (1. + x - x * << 
138   }                                            << 
139   else                                         << 
140   {                                            << 
141     res = G4Exp(G4Log(a) * onethird);          << 
142   }                                            << 
143   res = invert ? 1. / res : res;               << 
144   return res;                                  << 
145 }                                              << 
146                                                << 
147 // ------------------------------------------- << 
148                                                << 
149 G4double G4Pow::A13Low(const G4double a, const << 
150 {                                              << 
151   G4double res;                                << 
152   const auto  i    = static_cast<G4int>(4. * ( << 
153   const G4double y = 0.25 * i;                 << 
154   const G4double x = (a / y - 1.) * onethird;  << 
155   res              = lowa13[i] * (1. + x - x * << 
156   res              = invert ? 1. / res : res;  << 
157   return res;                                  << 
158 }                                              << 
159                                                   125 
160 // -------------------------------------------    126 // -------------------------------------------------------------------
161                                                   127 
162 G4double G4Pow::powN(G4double x, G4int n) cons    128 G4double G4Pow::powN(G4double x, G4int n) const
163 {                                                 129 {
164   if(0.0 == x)                                 << 130   if(0.0 == x)        { return 0.0; }
165   {                                            << 131   if(std::abs(n) > 8) { return std::pow(x, G4double(n)); }
166     return 0.0;                                << 
167   }                                            << 
168   if(std::abs(n) > 8)                          << 
169   {                                            << 
170     return std::pow(x, G4double(n));           << 
171   }                                            << 
172   G4double res = 1.0;                             132   G4double res = 1.0;
173   if(n >= 0)                                   << 133   if(n >= 0) { for(G4int i=0; i<n; ++i) { res *= x; } }
174   {                                            << 
175     for(G4int i = 0; i < n; ++i)               << 
176     {                                          << 
177       res *= x;                                << 
178     }                                          << 
179   }                                            << 
180   else if(n < 0)                                  134   else if(n < 0)
181   {                                               135   {
182     G4double y = 1.0 / x;                      << 136     G4double y = 1.0/x;
183     G4int nn   = -n;                           << 137     G4int nn = -n;
184     for(G4int i = 0; i < nn; ++i)              << 138     for(G4int i=0; i<nn; ++i) { res *= y; }
185     {                                          << 
186       res *= y;                                << 
187     }                                          << 
188   }                                               139   }
189   return res;                                     140   return res;
190 }                                                 141 }
191                                                   142