Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/adjoint/src/G4AdjointInterpolator.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 /processes/electromagnetic/adjoint/src/G4AdjointInterpolator.cc (Version 11.3.0) and /processes/electromagnetic/adjoint/src/G4AdjointInterpolator.cc (Version 11.0.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                                                    26 
 27 #include "G4AdjointInterpolator.hh"                27 #include "G4AdjointInterpolator.hh"
 28                                                    28 
 29 G4ThreadLocal G4AdjointInterpolator* G4Adjoint     29 G4ThreadLocal G4AdjointInterpolator* G4AdjointInterpolator::fInstance = nullptr;
 30                                                    30 
 31 //////////////////////////////////////////////     31 ///////////////////////////////////////////////////////
 32 G4AdjointInterpolator* G4AdjointInterpolator::     32 G4AdjointInterpolator* G4AdjointInterpolator::GetAdjointInterpolator()
 33 {                                                  33 {
 34   return GetInstance();                            34   return GetInstance();
 35 }                                                  35 }
 36                                                    36 
 37 //////////////////////////////////////////////     37 ///////////////////////////////////////////////////////
 38 G4AdjointInterpolator* G4AdjointInterpolator::     38 G4AdjointInterpolator* G4AdjointInterpolator::GetInstance()
 39 {                                                  39 {
 40   if(!fInstance)                                   40   if(!fInstance)
 41   {                                                41   {
 42     fInstance = new G4AdjointInterpolator;         42     fInstance = new G4AdjointInterpolator;
 43   }                                                43   }
 44   return fInstance;                                44   return fInstance;
 45 }                                                  45 }
 46                                                    46 
 47 //////////////////////////////////////////////     47 ///////////////////////////////////////////////////////
 48 G4AdjointInterpolator::G4AdjointInterpolator()     48 G4AdjointInterpolator::G4AdjointInterpolator() {}
 49                                                    49 
 50 //////////////////////////////////////////////     50 ///////////////////////////////////////////////////////
 51 G4AdjointInterpolator::~G4AdjointInterpolator(     51 G4AdjointInterpolator::~G4AdjointInterpolator() {}
 52                                                    52 
 53 //////////////////////////////////////////////     53 ///////////////////////////////////////////////////////
 54 G4double G4AdjointInterpolator::LinearInterpol     54 G4double G4AdjointInterpolator::LinearInterpolation(G4double& x, G4double& x1,
 55                                                    55                                                     G4double& x2, G4double& y1,
 56                                                    56                                                     G4double& y2)
 57 {                                                  57 {
 58   G4double res = y1 + (x - x1) * (y2 - y1) / (     58   G4double res = y1 + (x - x1) * (y2 - y1) / (x2 - x1);
 59   return res;                                      59   return res;
 60 }                                                  60 }
 61                                                    61 
 62 //////////////////////////////////////////////     62 ///////////////////////////////////////////////////////
 63 G4double G4AdjointInterpolator::LogarithmicInt     63 G4double G4AdjointInterpolator::LogarithmicInterpolation(
 64   G4double& x, G4double& x1, G4double& x2, G4d     64   G4double& x, G4double& x1, G4double& x2, G4double& y1, G4double& y2)
 65 {                                                  65 {
 66   if(y1 <= 0. || y2 <= 0. || x1 <= 0.)             66   if(y1 <= 0. || y2 <= 0. || x1 <= 0.)
 67     return LinearInterpolation(x, x1, x2, y1,      67     return LinearInterpolation(x, x1, x2, y1, y2);
 68   G4double B   = std::log(y2 / y1) / std::log(     68   G4double B   = std::log(y2 / y1) / std::log(x2 / x1);
 69   G4double A   = y1 / std::pow(x1, B);             69   G4double A   = y1 / std::pow(x1, B);
 70   G4double res = A * std::pow(x, B);               70   G4double res = A * std::pow(x, B);
 71   return res;                                      71   return res;
 72 }                                                  72 }
 73                                                    73 
 74 //////////////////////////////////////////////     74 ///////////////////////////////////////////////////////
 75 G4double G4AdjointInterpolator::ExponentialInt     75 G4double G4AdjointInterpolator::ExponentialInterpolation(
 76   G4double& x, G4double& x1, G4double& x2, G4d     76   G4double& x, G4double& x1, G4double& x2, G4double& y1, G4double& y2)
 77 {                                                  77 {
 78   G4double B   = (std::log(y2) - std::log(y1))     78   G4double B   = (std::log(y2) - std::log(y1)) / (x2 - x1);
 79   G4double A   = y1 * std::exp(-B * x1);           79   G4double A   = y1 * std::exp(-B * x1);
 80   G4double res = A * std::exp(B * x);              80   G4double res = A * std::exp(B * x);
 81   return res;                                      81   return res;
 82 }                                                  82 }
 83                                                    83 
 84 //////////////////////////////////////////////     84 ///////////////////////////////////////////////////////
 85 G4double G4AdjointInterpolator::Interpolation(     85 G4double G4AdjointInterpolator::Interpolation(G4double& x, G4double& x1,
 86                                                    86                                               G4double& x2, G4double& y1,
 87                                                    87                                               G4double& y2,
 88                                                <<  88                                               G4String InterPolMethod)
 89 {                                                  89 {
 90   if(InterPolMethod == "Log")                      90   if(InterPolMethod == "Log")
 91   {                                                91   {
 92     return LogarithmicInterpolation(x, x1, x2,     92     return LogarithmicInterpolation(x, x1, x2, y1, y2);
 93   }                                                93   }
 94   else if(InterPolMethod == "Lin")                 94   else if(InterPolMethod == "Lin")
 95   {                                                95   {
 96     return LinearInterpolation(x, x1, x2, y1,      96     return LinearInterpolation(x, x1, x2, y1, y2);
 97   }                                                97   }
 98   else if(InterPolMethod == "Exp")                 98   else if(InterPolMethod == "Exp")
 99   {                                                99   {
100     return ExponentialInterpolation(x, x1, x2,    100     return ExponentialInterpolation(x, x1, x2, y1, y2);
101   }                                               101   }
102   else                                            102   else
103   {                                               103   {
104     G4ExceptionDescription ed;                    104     G4ExceptionDescription ed;
105     ed << "The interpolation method that you i    105     ed << "The interpolation method that you invoked does not exist!\n";
106     G4Exception("G4AdjointInterpolator::Interp    106     G4Exception("G4AdjointInterpolator::Interpolation", "adoint001",
107                 FatalException, ed);              107                 FatalException, ed);
108     return 0.;                                    108     return 0.;
109   }                                               109   }
110 }                                                 110 }
111                                                   111 
112 //////////////////////////////////////////////    112 ///////////////////////////////////////////////////////
113 // only valid if x_vec is monotically increasi    113 // only valid if x_vec is monotically increasing
114 std::size_t G4AdjointInterpolator::FindPositio << 114 size_t G4AdjointInterpolator::FindPosition(G4double& x,
115                                            std << 115                                            std::vector<G4double>& x_vec, size_t,
116                                            std << 116                                            size_t)
117 {                                                 117 {
118   // most rapid method could be used probably     118   // most rapid method could be used probably
119                                                   119 
120   std::size_t ndim = x_vec.size();             << 120   size_t ndim = x_vec.size();
121   std::size_t ind1 = 0;                        << 121   size_t ind1 = 0;
122   std::size_t ind2 = ndim - 1;                 << 122   size_t ind2 = ndim - 1;
123                                                   123 
124   if(ndim > 1)                                    124   if(ndim > 1)
125   {                                               125   {
126     if(x_vec[0] < x_vec[1])                       126     if(x_vec[0] < x_vec[1])
127     {  // increasing                              127     {  // increasing
128       do                                          128       do
129       {                                           129       {
130         std::size_t midBin = (ind1 + ind2) / 2 << 130         size_t midBin = (ind1 + ind2) / 2;
131         if(x < x_vec[midBin])                     131         if(x < x_vec[midBin])
132           ind2 = midBin;                          132           ind2 = midBin;
133         else                                      133         else
134           ind1 = midBin;                          134           ind1 = midBin;
135       } while(ind2 - ind1 > 1);                   135       } while(ind2 - ind1 > 1);
136     }                                             136     }
137     else                                          137     else
138     {                                             138     {
139       do                                          139       do
140       {                                           140       {
141         std::size_t midBin = (ind1 + ind2) / 2 << 141         size_t midBin = (ind1 + ind2) / 2;
142         if(x < x_vec[midBin])                     142         if(x < x_vec[midBin])
143           ind1 = midBin;                          143           ind1 = midBin;
144         else                                      144         else
145           ind2 = midBin;                          145           ind2 = midBin;
146       } while(ind2 - ind1 > 1);                   146       } while(ind2 - ind1 > 1);
147     }                                             147     }
148   }                                               148   }
149                                                   149 
150   return ind1;                                    150   return ind1;
151 }                                                 151 }
152                                                   152 
153 //////////////////////////////////////////////    153 ///////////////////////////////////////////////////////
154 // only valid if x_vec is monotically increasi    154 // only valid if x_vec is monotically increasing
155 std::size_t G4AdjointInterpolator::FindPositio << 155 size_t G4AdjointInterpolator::FindPositionForLogVector(
156   G4double& log_x, std::vector<G4double>& log_    156   G4double& log_x, std::vector<G4double>& log_x_vec)
157 {                                                 157 {
158   // most rapid method could be used probably     158   // most rapid method could be used probably
159   return FindPosition(log_x, log_x_vec);          159   return FindPosition(log_x, log_x_vec);
160 }                                                 160 }
161                                                   161 
162 //////////////////////////////////////////////    162 ///////////////////////////////////////////////////////
163 G4double G4AdjointInterpolator::Interpolate(G4    163 G4double G4AdjointInterpolator::Interpolate(G4double& x,
164                                             st    164                                             std::vector<G4double>& x_vec,
165                                             st    165                                             std::vector<G4double>& y_vec,
166                                             co << 166                                             G4String InterPolMethod)
167 {                                                 167 {
168   std::size_t i = FindPosition(x, x_vec);      << 168   size_t i = FindPosition(x, x_vec);
169   return Interpolation(x, x_vec[i], x_vec[i +     169   return Interpolation(x, x_vec[i], x_vec[i + 1], y_vec[i], y_vec[i + 1],
170                        InterPolMethod);           170                        InterPolMethod);
171 }                                                 171 }
172                                                   172 
173 //////////////////////////////////////////////    173 ///////////////////////////////////////////////////////
174 G4double G4AdjointInterpolator::InterpolateWit    174 G4double G4AdjointInterpolator::InterpolateWithIndexVector(
175   G4double& x, std::vector<G4double>& x_vec, s    175   G4double& x, std::vector<G4double>& x_vec, std::vector<G4double>& y_vec,
176   std::vector<std::size_t>& index_vec, G4doubl << 176   std::vector<size_t>& index_vec, G4double x0,
177   G4double dx)  // only linear interpolation p    177   G4double dx)  // only linear interpolation possible
178 {                                                 178 {
179   std::size_t ind = 0;                         << 179   size_t ind = 0;
180   if(x > x0)                                      180   if(x > x0)
181     ind = G4int((x - x0) / dx);                << 181     ind = int((x - x0) / dx);
182   if(ind >= index_vec.size() - 1)                 182   if(ind >= index_vec.size() - 1)
183     ind = index_vec.size() - 2;                   183     ind = index_vec.size() - 2;
184   std::size_t ind1 = index_vec[ind];           << 184   size_t ind1 = index_vec[ind];
185   std::size_t ind2 = index_vec[ind + 1];       << 185   size_t ind2 = index_vec[ind + 1];
186   if(ind1 > ind2)                                 186   if(ind1 > ind2)
187   {                                               187   {
188     std::size_t ind11 = ind1;                  << 188     size_t ind11 = ind1;
189     ind1         = ind2;                          189     ind1         = ind2;
190     ind2         = ind11;                         190     ind2         = ind11;
191   }                                               191   }
192   ind = FindPosition(x, x_vec, ind1, ind2);       192   ind = FindPosition(x, x_vec, ind1, ind2);
193   return Interpolation(x, x_vec[ind], x_vec[in    193   return Interpolation(x, x_vec[ind], x_vec[ind + 1], y_vec[ind],
194                        y_vec[ind + 1], "Lin");    194                        y_vec[ind + 1], "Lin");
195 }                                                 195 }
196                                                   196 
197 //////////////////////////////////////////////    197 ///////////////////////////////////////////////////////
198 G4double G4AdjointInterpolator::InterpolateFor    198 G4double G4AdjointInterpolator::InterpolateForLogVector(
199   G4double& log_x, std::vector<G4double>& log_    199   G4double& log_x, std::vector<G4double>& log_x_vec,
200   std::vector<G4double>& log_y_vec)               200   std::vector<G4double>& log_y_vec)
201 {                                                 201 {
202   std::size_t i = FindPositionForLogVector(log << 202   size_t i = FindPositionForLogVector(log_x, log_x_vec);
203                                                   203 
204   G4double log_y = LinearInterpolation(log_x,     204   G4double log_y = LinearInterpolation(log_x, log_x_vec[i], log_x_vec[i + 1],
205                                        log_y_v    205                                        log_y_vec[i], log_y_vec[i + 1]);
206   return log_y;                                   206   return log_y;
207 }                                                 207 }
208                                                   208