include/Geant4/G4PhysicsVector.icc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // G4PhysicsVector inline methods implementation
0027 //
0028 // Authors:
0029 // - 02 Dec. 1995, G.Cosmo: Structure created based on object model
0030 // - 03 Mar. 1996, K.Amako: Implemented the 1st version
0031 // --------------------------------------------------------------------
0032 inline G4double G4PhysicsVector::operator[](const std::size_t index) const
0033 {
0034   return dataVector[index];
0035 }
0036
0037 // ---------------------------------------------------------------
0038 inline G4double G4PhysicsVector::operator()(const std::size_t index) const
0039 {
0040   return dataVector[index];
0041 }
0042
0043 // ---------------------------------------------------------------
0044 inline G4double G4PhysicsVector::Energy(const std::size_t index) const
0045 {
0046   return binVector[index];
0047 }
0048
0049 // ---------------------------------------------------------------
0050 inline G4double
0051 G4PhysicsVector::GetLowEdgeEnergy(const std::size_t index) const
0052 {
0053   return binVector[index];
0054 }
0055
0056 // ---------------------------------------------------------------
0057 inline G4double G4PhysicsVector::GetMinEnergy() const
0058 {
0059   return edgeMin;
0060 }
0061
0062 // ---------------------------------------------------------------
0063 inline G4double G4PhysicsVector::GetMaxEnergy() const
0064 {
0065   return edgeMax;
0066 }
0067
0068 // ---------------------------------------------------------------
0069 inline G4double G4PhysicsVector::GetMinValue() const
0070 {
0071   return (numberOfNodes > 0) ? dataVector[0] : 0.0;
0072 }
0073
0074 // ---------------------------------------------------------------
0075 inline G4double G4PhysicsVector::GetMaxValue() const
0076 {
0077   return (numberOfNodes > 0) ? dataVector[numberOfNodes - 1] : 0.0;
0078 }
0079
0080 // ---------------------------------------------------------------
0081 inline std::size_t G4PhysicsVector::GetVectorLength() const
0082 {
0083   return numberOfNodes;
0084 }
0085
0086 // ---------------------------------------------------------------
0087 inline void G4PhysicsVector::PutValue(std::size_t index, G4double theValue)
0088 {
0089   if(index >= numberOfNodes)
0090   {
0091     PrintPutValueError(index, theValue, "PutValue(..) ");
0092   }
0093   else
0094   {
0095     dataVector[index] = theValue;
0096   }
0097 }
0098
0099 // ---------------------------------------------------------------
0100 inline G4PhysicsVectorType G4PhysicsVector::GetType() const
0101 {
0102   return type;
0103 }
0104
0105 // ---------------------------------------------------------------
0106 inline G4bool G4PhysicsVector::GetSpline() const
0107 {
0108   return useSpline;
0109 }
0110
0111 // ---------------------------------------------------------------
0112 inline void G4PhysicsVector::SetVerboseLevel(G4int value)
0113 {
0114   verboseLevel = value;
0115 }
0116
0117 // ---------------------------------------------------------------
0118 inline G4double
0119 G4PhysicsVector::FindLinearEnergy(const G4double rand) const
0120 {
0121   return GetEnergy(rand*dataVector[numberOfNodes - 1]);
0122 }
0123
0124 // ---------------------------------------------------------------
0125 inline G4double G4PhysicsVector::Interpolation(const std::size_t idx,
0126                                                const G4double e) const
0127 {
0128   // perform the interpolation
0129   const G4double x1 = binVector[idx];
0130   const G4double dl = binVector[idx + 1] - x1;
0131
0132   const G4double y1 = dataVector[idx];
0133   const G4double dy = dataVector[idx + 1] - y1;
0134
0135   // note: all corner cases of the previous methods are covered and eventually
0136   //       gives b=0/1 that results in y=y0\y_{N-1} if e<=x[0]/e>=x[N-1] or
0137   //       y=y_i/y_{i+1} if e<x[i]/e>=x[i+1] due to small numerical errors
0138   const G4double b = (e - x1) / dl;
0139
0140   G4double res = y1 + b * dy;
0141
0142   if (useSpline)  // spline interpolation
0143   {
0144     const G4double c0 = (2.0 - b) * secDerivative[idx];
0145     const G4double c1 = (1.0 + b) * secDerivative[idx + 1];
0146     res += (b * (b - 1.0)) * (c0 + c1) * (dl * dl * (1.0/6.0));
0147   }
0148
0149   return res;
0150 }
0151
0152 // ---------------------------------------------------------------
0153 inline std::size_t G4PhysicsVector::ComputeLogVectorBin(
0154   const G4double loge) const
0155 {
0156   return static_cast<std::size_t>( std::min( static_cast<G4int>((loge - logemin) * invdBin),
0157     static_cast<G4int>(idxmax) ) );
0158 }
0159
0160 // ---------------------------------------------------------------
0161 inline std::size_t
0162 G4PhysicsVector::LogBin(const G4double e, const G4double loge) const
0163 {
0164   std::size_t idx =
0165   scale[std::min( static_cast<G4int>((loge - lmin1) * iBin1),
0166                   static_cast<G4int>(imax1) )];
0167   for (; idx <= idxmax; ++idx)
0168   {
0169     if (e >= binVector[idx] && e <= binVector[idx + 1]) { break; }
0170   }
0171   return idx;
0172 }
0173
0174 // ---------------------------------------------------------------
0175 inline std::size_t G4PhysicsVector::BinaryBin(const G4double e) const
0176 {
0177   // Bin location proposed by K.Genser (FNAL)
0178   return std::lower_bound(binVector.cbegin(), binVector.cend(), e) -
0179     binVector.cbegin() - 1;
0180 }
0181
0182 // ---------------------------------------------------------------
0183 inline std::size_t G4PhysicsVector::GetBin(const G4double e) const
0184 {
0185   std::size_t bin;
0186   switch(type)
0187   {
0188     case T_G4PhysicsLogVector:
0189       bin = ComputeLogVectorBin(G4Log(e));
0190       break;
0191
0192     case T_G4PhysicsLinearVector:
0193       bin = static_cast<std::size_t>( std::min( static_cast<G4int>((e - edgeMin) * invdBin),
0194                                                 static_cast<G4int>(idxmax) ) );
0195       break;
0196
0197     default:
0198       bin = (nLogNodes > 0) ? LogBin(e, G4Log(e)) : BinaryBin(e);
0199   }
0200   return bin;
0201 }
0202
0203 // ---------------------------------------------------------------
0204 inline G4double
0205 G4PhysicsVector::Value(const G4double e, std::size_t& idx) const
0206 {
0207   G4double res;
0208   if (idx + 1 < numberOfNodes &&
0209       e >= binVector[idx] && e <= binVector[idx+1])
0210   {
0211     res = Interpolation(idx, e);
0212   }
0213   else if (e > edgeMin && e < edgeMax)
0214   {
0215     idx = GetBin(e);
0216     res = Interpolation(idx, e);
0217   }
0218   else if(e <= edgeMin)
0219   {
0220     res = dataVector[0];
0221     idx = 0;
0222   }
0223   else
0224   {
0225     res = dataVector[idxmax + 1];
0226     idx = idxmax;
0227   }
0228   return res;
0229 }
0230
0231 // ---------------------------------------------------------------
0232 inline G4double G4PhysicsVector::Value(G4double e) const
0233 {
0234   G4double res;
0235   if (e > edgeMin && e < edgeMax)
0236   {
0237     const std::size_t idx = GetBin(e);
0238     res = Interpolation(idx, e);
0239   }
0240   else if(e <= edgeMin)
0241   {
0242     res = dataVector[0];
0243   }
0244   else
0245   {
0246     res = dataVector[idxmax + 1];
0247   }
0248   return res;
0249 }
0250
0251 // ---------------------------------------------------------------
0252 inline G4double G4PhysicsVector::GetValue(G4double e, G4bool&) const
0253 {
0254   return Value(e);
0255 }
0256
0257 // ---------------------------------------------------------------
0258 inline G4double
0259 G4PhysicsVector::LogVectorValue(const G4double e, const G4double loge) const
0260 {
0261   G4double res;
0262   if (e > edgeMin && e < edgeMax)
0263   {
0264     const std::size_t idx = ComputeLogVectorBin(loge);
0265     res = Interpolation(idx, e);
0266   }
0267   else if (e <= edgeMin)
0268   {
0269     res = dataVector[0];
0270   }
0271   else
0272   {
0273     res = dataVector[idxmax - 1];
0274   }
0275   return res;
0276 }
0277
0278 // ---------------------------------------------------------------
0279 inline G4double
0280 G4PhysicsVector::LogFreeVectorValue(const G4double e, const G4double loge) const
0281 {
0282   G4double res;
0283   if (e > edgeMin && e < edgeMax)
0284   {
0285     const std::size_t idx = LogBin(e, loge);
0286     res = Interpolation(idx, e);
0287   }
0288   else if (e <= edgeMin)
0289   {
0290     res = dataVector[0];
0291   }
0292   else
0293   {
0294     res = dataVector[idxmax + 1];
0295   }
0296   return res;
0297 }
0298
0299 // ---------------------------------------------------------------