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 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 class implemenation                   <<  26 // $Id$
 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 //
                                                   >>  42 // -------------------------------------------------------------------
 33                                                    43 
 34 #include "G4Pow.hh"                                44 #include "G4Pow.hh"
 35 #include "G4Threading.hh"                      << 
 36                                                    45 
 37 G4Pow* G4Pow::fpInstance = nullptr;            <<  46 G4Pow* G4Pow::fInstance = 0;
 38                                                    47 
 39 // -------------------------------------------     48 // -------------------------------------------------------------------
 40                                                    49 
 41 G4Pow* G4Pow::GetInstance()                        50 G4Pow* G4Pow::GetInstance()
 42 {                                                  51 {
 43   if(fpInstance == nullptr)                    <<  52   if (fInstance == 0)
 44   {                                                53   {
 45     static G4Pow geant4pow;                    <<  54     static G4Pow manager;
 46     fpInstance = &geant4pow;                   <<  55     fInstance = &manager;
 47   }                                                56   }
 48   return fpInstance;                           <<  57   return fInstance;
 49 }                                                  58 }
 50                                                    59 
 51 // -------------------------------------------     60 // -------------------------------------------------------------------
 52                                                    61 
 53 G4Pow::G4Pow()                                     62 G4Pow::G4Pow()
 54 {                                                  63 {
 55 #ifdef G4MULTITHREADED                         <<  64   const G4int maxZ = 512; 
 56   if(G4Threading::IsWorkerThread())            <<  65   const G4int maxZfact = 170; 
 57   {                                            << 
 58     G4Exception("G4Pow::G4Pow()", "InvalidSetu << 
 59                 "Attempt to instantiate G4Pow  << 
 60   }                                            << 
 61 #endif                                         << 
 62   const G4int maxZ     = 512;                  << 
 63   const G4int maxZfact = 170;                  << 
 64   const G4int numLowA  = 17;                   << 
 65                                                << 
 66   maxA    = -0.6 + maxZ;                       << 
 67   maxLowA = 4.0;                               << 
 68   maxA2   = 1.25 + max2 * 0.2;                 << 
 69   maxAexp = -0.76 + maxZfact * 0.5;            << 
 70                                                << 
 71   ener.resize(max2 + 1, 1.0);                  << 
 72   logen.resize(max2 + 1, 0.0);                 << 
 73   lz2.resize(max2 + 1, 0.0);                   << 
 74   pz13.resize(maxZ, 0.0);                      << 
 75   lowa13.resize(numLowA, 0.0);                 << 
 76   lz.resize(maxZ, 0.0);                        << 
 77   fexp.resize(maxZfact, 0.0);                  << 
 78   fact.resize(maxZfact, 0.0);                  << 
 79   logfact.resize(maxZ, 0.0);                   << 
 80                                                    66 
 81   G4double f    = 1.0;                         <<  67   pz13.resize(maxZ,0.0);
 82   G4double logf = 0.0;                         <<  68   lz.resize(maxZ,0.0);
 83   fact[0]       = 1.0;                         <<  69   fact.resize(maxZfact,0.0);
 84   fexp[0]       = 1.0;                         <<  70   logfact.resize(maxZ,0.0);
 85                                                    71 
 86   for(G4int i = 1; i <= max2; ++i)             <<  72   onethird = 1.0/3.0;
 87   {                                            <<  73   G4double f = 1.0;
 88     ener[i]  = powN(500., i);                  <<  74   G4double logf = 0.0;
 89     logen[i] = G4Log(ener[i]);                 <<  75   fact[0] = 1.0;
 90     lz2[i]   = G4Log(1.0 + i * 0.2);           << 
 91   }                                            << 
 92                                                    76 
 93   for(G4int i = 1; i < maxZ; ++i)              <<  77   for(G4int i=1; i<maxZ; ++i)
 94   {                                                78   {
 95     auto  x = G4double(i);                     <<  79     G4double x  = G4double(i);
 96     pz13[i]    = std::pow(x, onethird);        <<  80     pz13[i] = std::pow(x,onethird);
 97     lz[i]      = G4Log(x);                     <<  81     lz[i]   = std::log(x);
 98     if(i < maxZfact)                           <<  82     if(i < maxZfact) { 
 99     {                                          <<  83       f *= x; 
100       f *= x;                                  << 
101       fact[i] = f;                                 84       fact[i] = f;
102       fexp[i] = G4Exp(0.5 * i);                << 
103     }                                              85     }
104     logf += lz[i];                                 86     logf += lz[i];
105     logfact[i] = logf;                             87     logfact[i] = logf;
106   }                                                88   }
107                                                << 
108   for(G4int i = 4; i < numLowA; ++i)           << 
109   {                                            << 
110     lowa13[i] = std::pow(0.25 * i, onethird);  << 
111   }                                            << 
112 }                                              << 
113                                                << 
114 // ------------------------------------------- << 
115                                                << 
116 G4double G4Pow::A13(G4double A) const          << 
117 {                                              << 
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 }                                                  89 }
146                                                    90 
147 // -------------------------------------------     91 // -------------------------------------------------------------------
148                                                    92 
149 G4double G4Pow::A13Low(const G4double a, const <<  93 G4Pow::~G4Pow()
150 {                                              <<  94 {}
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                                                    95 
160 // -------------------------------------------     96 // -------------------------------------------------------------------
161                                                    97 
162 G4double G4Pow::powN(G4double x, G4int n) cons <<  98 G4double G4Pow::powN(G4double x, G4int n)
163 {                                                  99 {
164   if(0.0 == x)                                 << 100   if(std::abs(n) > 8) { return std::pow(x, G4double(n)); }
165   {                                            << 
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;                             101   G4double res = 1.0;
173   if(n >= 0)                                   << 102   if(n >= 0) { for(G4int i=0; i<n; ++i) { res *= x; } }
                                                   >> 103   else if((n < 0) && (x != 0.0))
174   {                                               104   {
175     for(G4int i = 0; i < n; ++i)               << 105     G4double y = 1.0/x;
176     {                                          << 106     G4int nn = -n;
177       res *= x;                                << 107     for(G4int i=0; i<nn; ++i) { res *= y; }
178     }                                          << 
179   }                                            << 
180   else if(n < 0)                               << 
181   {                                            << 
182     G4double y = 1.0 / x;                      << 
183     G4int nn   = -n;                           << 
184     for(G4int i = 0; i < nn; ++i)              << 
185     {                                          << 
186       res *= y;                                << 
187     }                                          << 
188   }                                               108   }
189   return res;                                     109   return res;
190 }                                                 110 }
191                                                   111