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.0.p3)


  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 74256 2013-10-02 14:24:02Z 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"
 35 #include "G4Threading.hh"                          47 #include "G4Threading.hh"
 36                                                    48 
 37 G4Pow* G4Pow::fpInstance = nullptr;            <<  49 G4Pow* G4Pow::fpInstance = 0;
 38                                                    50 
 39 // -------------------------------------------     51 // -------------------------------------------------------------------
 40                                                    52 
 41 G4Pow* G4Pow::GetInstance()                        53 G4Pow* G4Pow::GetInstance()
 42 {                                                  54 {
 43   if(fpInstance == nullptr)                    <<  55   if (fpInstance == 0)
 44   {                                                56   {
 45     static G4Pow geant4pow;                    <<  57     fpInstance = new G4Pow;
 46     fpInstance = &geant4pow;                   << 
 47   }                                                58   }
 48   return fpInstance;                               59   return fpInstance;
 49 }                                                  60 }
 50                                                    61 
 51 // -------------------------------------------     62 // -------------------------------------------------------------------
 52                                                    63 
 53 G4Pow::G4Pow()                                     64 G4Pow::G4Pow()
                                                   >>  65   : onethird(1.0/3.0), max2(5)
 54 {                                                  66 {
 55 #ifdef G4MULTITHREADED                         <<  67   if(G4Threading::IsWorkerThread() == true) { 
 56   if(G4Threading::IsWorkerThread())            <<  68     G4Exception ("G4Pow::G4Pow()", "glob090", FatalException, 
 57   {                                            <<  69                  "Attempt to instantiate G4Pow in worker thread");
 58     G4Exception("G4Pow::G4Pow()", "InvalidSetu << 
 59                 "Attempt to instantiate G4Pow  << 
 60   }                                                70   }
 61 #endif                                         <<  71 
 62   const G4int maxZ     = 512;                  <<  72   const G4int maxZ = 512; 
 63   const G4int maxZfact = 170;                  <<  73   const G4int maxZfact = 170; 
 64   const G4int numLowA  = 17;                   << 
 65                                                    74 
 66   maxA    = -0.6 + maxZ;                           75   maxA    = -0.6 + maxZ;
 67   maxLowA = 4.0;                               <<  76   maxA2   = 1.25 + max2*0.2;
 68   maxA2   = 1.25 + max2 * 0.2;                 <<  77   maxAexp = -0.76+ maxZfact*0.5;
 69   maxAexp = -0.76 + maxZfact * 0.5;            <<  78 
 70                                                <<  79   ener.resize(max2+1,1.0);
 71   ener.resize(max2 + 1, 1.0);                  <<  80   logen.resize(max2+1,0.0);
 72   logen.resize(max2 + 1, 0.0);                 <<  81   lz2.resize(max2+1,0.0);
 73   lz2.resize(max2 + 1, 0.0);                   <<  82   pz13.resize(maxZ,0.0);
 74   pz13.resize(maxZ, 0.0);                      <<  83   lz.resize(maxZ,0.0);
 75   lowa13.resize(numLowA, 0.0);                 <<  84   fexp.resize(maxZfact,0.0);
 76   lz.resize(maxZ, 0.0);                        <<  85   fact.resize(maxZfact,0.0);
 77   fexp.resize(maxZfact, 0.0);                  <<  86   logfact.resize(maxZ,0.0);
 78   fact.resize(maxZfact, 0.0);                  << 
 79   logfact.resize(maxZ, 0.0);                   << 
 80                                                    87 
 81   G4double f    = 1.0;                         <<  88   G4double f = 1.0;
 82   G4double logf = 0.0;                             89   G4double logf = 0.0;
 83   fact[0]       = 1.0;                         <<  90   fact[0] = 1.0;
 84   fexp[0]       = 1.0;                         <<  91   fexp[0] = 1.0;
 85                                                    92 
 86   for(G4int i = 1; i <= max2; ++i)             <<  93   for(G4int i=1; i<=max2; ++i)
 87   {                                                94   {
 88     ener[i]  = powN(500., i);                  <<  95     ener[i] = powN(500.,i); 
 89     logen[i] = G4Log(ener[i]);                 <<  96     logen[i]= G4Log(ener[i]); 
 90     lz2[i]   = G4Log(1.0 + i * 0.2);           <<  97     lz2[i]  = G4Log(1.0 + i*0.2);
 91   }                                                98   }
 92                                                    99 
 93   for(G4int i = 1; i < maxZ; ++i)              << 100   for(G4int i=1; i<maxZ; ++i)
 94   {                                               101   {
 95     auto  x = G4double(i);                     << 102     G4double x  = G4double(i);
 96     pz13[i]    = std::pow(x, onethird);        << 103     pz13[i] = std::pow(x,onethird);
 97     lz[i]      = G4Log(x);                     << 104     lz[i]   = G4Log(x);
 98     if(i < maxZfact)                              105     if(i < maxZfact)
 99     {                                          << 106     { 
100       f *= x;                                  << 107       f *= x; 
101       fact[i] = f;                                108       fact[i] = f;
102       fexp[i] = G4Exp(0.5 * i);                << 109       fexp[i] = G4Exp(0.5*i);
103     }                                             110     }
104     logf += lz[i];                                111     logf += lz[i];
105     logfact[i] = logf;                            112     logfact[i] = logf;
106   }                                               113   }
107                                                << 
108   for(G4int i = 4; i < numLowA; ++i)           << 
109   {                                            << 
110     lowa13[i] = std::pow(0.25 * i, onethird);  << 
111   }                                            << 
112 }                                                 114 }
113                                                   115 
114 // -------------------------------------------    116 // -------------------------------------------------------------------
115                                                   117 
116 G4double G4Pow::A13(G4double A) const          << 118 G4Pow::~G4Pow()
117 {                                                 119 {
118   G4double res = 0.;                           << 120   delete fpInstance; fpInstance = 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 }                                                 121 }
159                                                   122 
160 // -------------------------------------------    123 // -------------------------------------------------------------------
161                                                   124 
162 G4double G4Pow::powN(G4double x, G4int n) cons    125 G4double G4Pow::powN(G4double x, G4int n) const
163 {                                                 126 {
164   if(0.0 == x)                                 << 127   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;                             128   G4double res = 1.0;
173   if(n >= 0)                                   << 129   if(n >= 0) { for(G4int i=0; i<n; ++i) { res *= x; } }
                                                   >> 130   else if((n < 0) && (x != 0.0))
174   {                                               131   {
175     for(G4int i = 0; i < n; ++i)               << 132     G4double y = 1.0/x;
176     {                                          << 133     G4int nn = -n;
177       res *= x;                                << 134     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   }                                               135   }
189   return res;                                     136   return res;
190 }                                                 137 }
191                                                   138