include/Geant4/G4ParticleHPInterpolator.hh

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 //
0027 // 080809 Change interpolation scheme of "histogram", now using LinearLinear
0028 //        For multidimensional interpolations By T. Koi
0029 //
0030 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
0031 //
0032 #ifndef G4ParticleHPInterpolator_h
0033 #define G4ParticleHPInterpolator_h 1
0034
0035 #include "G4Exp.hh"
0036 #include "G4HadronicException.hh"
0037 #include "G4InterpolationScheme.hh"
0038 #include "G4Log.hh"
0039 #include "G4ios.hh"
0040 #include "Randomize.hh"
0041 #include "globals.hh"
0042
0043 class G4ParticleHPInterpolator
0044 {
0045   public:
0046     G4ParticleHPInterpolator() = default;
0047     ~G4ParticleHPInterpolator() = default;
0048
0049     inline G4double Lin(G4double x, G4double x1, G4double x2, G4double y1, G4double y2)
0050     {
0051       G4double slope = 0, off = 0;
0052       if (x2 - x1 == 0) return (y2 + y1) / 2.;
0053       slope = (y2 - y1) / (x2 - x1);
0054       off = y2 - x2 * slope;
0055       G4double y = x * slope + off;
0056       return y;
0057     }
0058
0059     inline G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2,
0060                                 G4double y1, G4double y2) const;
0061     inline G4double Interpolate2(G4InterpolationScheme aScheme, G4double x, G4double x1,
0062                                  G4double x2, G4double y1, G4double y2) const;
0063
0064     G4double GetBinIntegral(const G4InterpolationScheme& aScheme, const G4double x1,
0065                             const G4double x2, const G4double y1, const G4double y2);
0066
0067     G4double GetWeightedBinIntegral(const G4InterpolationScheme& aScheme, const G4double x1,
0068                                     const G4double x2, const G4double y1, const G4double y2);
0069
0070   private:
0071     inline G4double Histogram(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const;
0072     inline G4double LinearLinear(G4double x, G4double x1, G4double x2, G4double y1,
0073                                  G4double y2) const;
0074     inline G4double LinearLogarithmic(G4double x, G4double x1, G4double x2, G4double y1,
0075                                       G4double y2) const;
0076     inline G4double LogarithmicLinear(G4double x, G4double x1, G4double x2, G4double y1,
0077                                       G4double y2) const;
0078     inline G4double LogarithmicLogarithmic(G4double x, G4double x1, G4double x2, G4double y1,
0079                                            G4double y2) const;
0080     inline G4double Random(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const;
0081 };
0082
0083 inline G4double G4ParticleHPInterpolator::Interpolate(G4InterpolationScheme aScheme, G4double x,
0084                                                       G4double x1, G4double x2, G4double y1,
0085                                                       G4double y2) const
0086 {
0087   G4double result(0);
0088   G4int theScheme = aScheme;
0089   theScheme = theScheme % CSTART_;
0090   switch (theScheme) {
0091     case 1:
0092       // 080809
0093       // result = Histogram(x, x1, x2, y1, y2);
0094       result = LinearLinear(x, x1, x2, y1, y2);
0095       break;
0096     case 2:
0097       result = LinearLinear(x, x1, x2, y1, y2);
0098       break;
0099     case 3:
0100       result = LinearLogarithmic(x, x1, x2, y1, y2);
0101       break;
0102     case 4:
0103       result = LogarithmicLinear(x, x1, x2, y1, y2);
0104       break;
0105     case 5:
0106       result = LogarithmicLogarithmic(x, x1, x2, y1, y2);
0107       break;
0108     case 6:
0109       result = Random(x, x1, x2, y1, y2);
0110       break;
0111     default:
0112       G4cout << "theScheme = " << theScheme << G4endl;
0113       throw G4HadronicException(__FILE__, __LINE__,
0114                                 "G4ParticleHPInterpolator::Carthesian Invalid InterpolationScheme");
0115       break;
0116   }
0117   return result;
0118 }
0119
0120 inline G4double G4ParticleHPInterpolator::Interpolate2(G4InterpolationScheme aScheme, G4double x,
0121                                                        G4double x1, G4double x2, G4double y1,
0122                                                        G4double y2) const
0123 {
0124   G4double result(0);
0125   G4int theScheme = aScheme;
0126   theScheme = theScheme % CSTART_;
0127   switch (theScheme) {
0128     case 1:
0129       result = Histogram(x, x1, x2, y1, y2);
0130       break;
0131     case 2:
0132       result = LinearLinear(x, x1, x2, y1, y2);
0133       break;
0134     case 3:
0135       result = LinearLogarithmic(x, x1, x2, y1, y2);
0136       break;
0137     case 4:
0138       result = LogarithmicLinear(x, x1, x2, y1, y2);
0139       break;
0140     case 5:
0141       result = LogarithmicLogarithmic(x, x1, x2, y1, y2);
0142       break;
0143     case 6:
0144       result = Random(x, x1, x2, y1, y2);
0145       break;
0146     default:
0147       G4cout << "theScheme = " << theScheme << G4endl;
0148       throw G4HadronicException(__FILE__, __LINE__,
0149                                 "G4ParticleHPInterpolator::Carthesian Invalid InterpolationScheme");
0150       break;
0151   }
0152   return result;
0153 }
0154
0155 inline G4double G4ParticleHPInterpolator::Histogram(G4double, G4double, G4double, G4double y1,
0156                                                     G4double) const
0157 {
0158   G4double result;
0159   result = y1;
0160   return result;
0161 }
0162
0163 inline G4double G4ParticleHPInterpolator::LinearLinear(G4double x, G4double x1, G4double x2,
0164                                                        G4double y1, G4double y2) const
0165 {
0166   G4double slope = 0, off = 0;
0167   if (x2 - x1 == 0) return (y2 + y1) / 2.;
0168   slope = (y2 - y1) / (x2 - x1);
0169   off = y2 - x2 * slope;
0170   G4double y = x * slope + off;
0171   return y;
0172 }
0173
0174 inline G4double G4ParticleHPInterpolator::LinearLogarithmic(G4double x, G4double x1, G4double x2,
0175                                                             G4double y1, G4double y2) const
0176 {
0177   G4double result;
0178   if (x == 0)
0179     result = y1 + y2 / 2.;
0180   else if (x1 == 0)
0181     result = y1;
0182   else if (x2 == 0)
0183     result = y2;
0184   else
0185     result = LinearLinear(G4Log(x), G4Log(x1), G4Log(x2), y1, y2);
0186   return result;
0187 }
0188
0189 inline G4double G4ParticleHPInterpolator::LogarithmicLinear(G4double x, G4double x1, G4double x2,
0190                                                             G4double y1, G4double y2) const
0191 {
0192   G4double result;
0193   if (y1 == 0 || y2 == 0)
0194     result = 0;
0195   else {
0196     result = LinearLinear(x, x1, x2, G4Log(y1), G4Log(y2));
0197     result = G4Exp(result);
0198   }
0199   return result;
0200 }
0201
0202 inline G4double G4ParticleHPInterpolator::LogarithmicLogarithmic(G4double x, G4double x1,
0203                                                                  G4double x2, G4double y1,
0204                                                                  G4double y2) const
0205 {
0206   if (x == 0) return y1 + y2 / 2.;
0207   if (x1 == 0) return y1;
0208   if (x2 == 0) return y2;
0209   G4double result;
0210   if (y1 == 0 || y2 == 0)
0211     result = 0;
0212   else {
0213     result = LinearLinear(G4Log(x), G4Log(x1), G4Log(x2), G4Log(y1), G4Log(y2));
0214     result = G4Exp(result);
0215   }
0216   return result;
0217 }
0218
0219 inline G4double G4ParticleHPInterpolator::Random(G4double, G4double, G4double, G4double y1,
0220                                                  G4double y2) const
0221 {
0222   G4double result;
0223   result = y1 + G4UniformRand() * (y2 - y1);
0224   return result;
0225 }
0226
0227 #endif