Geant4_10
G4NeutronHPInterpolator.hh
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 //
28 // 080809 Change interpolation scheme of "histogram", now using LinearLinear
29 // For multidimensional interpolations By T. Koi
30 //
31 #ifndef G4NeutronHPInterpolator_h
32 #define G4NeutronHPInterpolator_h 1
33 
34 #include "globals.hh"
35 #include "G4InterpolationScheme.hh"
36 #include "Randomize.hh"
37 #include "G4ios.hh"
38 #include "G4HadronicException.hh"
39 
40 
42 {
43  public:
44 
47  {
48  // G4cout <<"deleted the interpolator"<<G4endl;
49  }
50 
52  {
53  G4double slope=0, off=0;
54  if(x2-x1==0) return (y2+y1)/2.;
55  slope = (y2-y1)/(x2-x1);
56  off = y2-x2*slope;
57  G4double y = x*slope+off;
58  return y;
59  }
60 
63  G4double y1, G4double y2) const;
64 
65  G4double
66  GetBinIntegral(const G4InterpolationScheme & aScheme,
67  const G4double x1,const G4double x2,const G4double y1,const G4double y2);
68 
69  G4double
71  const G4double x1,const G4double x2,const G4double y1,const G4double y2);
72 
73  private:
74 
75  inline G4double Histogram(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const;
76  inline G4double LinearLinear(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const;
77  inline G4double LinearLogarithmic(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const;
78  inline G4double LogarithmicLinear(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const;
79  inline G4double LogarithmicLogarithmic(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const;
80  inline G4double Random(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const;
81 
82 };
83 
87 {
88  G4double result(0);
89  G4int theScheme = aScheme;
90  theScheme = theScheme%CSTART_;
91  switch(theScheme)
92  {
93  case 1:
94  //080809
95  //result = Histogram(x, x1, x2, y1, y2);
96  result = LinearLinear(x, x1, x2, y1, y2);
97  break;
98  case 2:
99  result = LinearLinear(x, x1, x2, y1, y2);
100  break;
101  case 3:
102  result = LinearLogarithmic(x, x1, x2, y1, y2);
103  break;
104  case 4:
105  result = LogarithmicLinear(x, x1, x2, y1, y2);
106  break;
107  case 5:
108  result = LogarithmicLogarithmic(x, x1, x2, y1, y2);
109  break;
110  case 6:
111  result = Random(x, x1, x2, y1, y2);
112  break;
113  default:
114  G4cout << "theScheme = "<<theScheme<<G4endl;
115  throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPInterpolator::Carthesian Invalid InterpolationScheme");
116  break;
117  }
118  return result;
119 }
120 
121 inline G4double G4NeutronHPInterpolator::
122 Histogram(G4double , G4double , G4double , G4double y1, G4double ) const
123 {
125  result = y1;
126  return result;
127 }
128 
129 inline G4double G4NeutronHPInterpolator::
130 LinearLinear(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
131 {
132  G4double slope=0, off=0;
133  if(x2-x1==0) return (y2+y1)/2.;
134  slope = (y2-y1)/(x2-x1);
135  off = y2-x2*slope;
136  G4double y = x*slope+off;
137  return y;
138 }
139 
140 inline G4double G4NeutronHPInterpolator::
141 LinearLogarithmic(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
142 {
144  if(x==0) result = y1+y2/2.;
145  else if(x1==0) result = y1;
146  else if(x2==0) result = y2;
147  else result = LinearLinear(std::log(x), std::log(x1), std::log(x2), y1, y2);
148  return result;
149 }
150 
151 inline G4double G4NeutronHPInterpolator::
152 LogarithmicLinear(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
153 {
155  if(y1==0||y2==0) result = 0;
156  else
157  {
158  result = LinearLinear(x, x1, x2, std::log(y1), std::log(y2));
159  result = std::exp(result);
160  }
161  return result;
162 }
163 
164 inline G4double G4NeutronHPInterpolator::
165 LogarithmicLogarithmic(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
166 {
168  if(x==0) result = y1+y2/2.;
169  else if(x1==0) result = y1;
170  else if(x2==0) result = y2;
171  if(y1==0||y2==0) result = 0;
172  else
173  {
174  result = LinearLinear(std::log(x), std::log(x1), std::log(x2), std::log(y1), std::log(y2));
175  result = std::exp(result);
176  }
177  return result;
178 }
179 
180 inline G4double G4NeutronHPInterpolator::
181 Random(G4double , G4double , G4double , G4double y1, G4double y2) const
182 {
184  result = y1+G4UniformRand()*(y2-y1);
185  return result;
186 }
187 
188 #endif
Double_t y2[nxs]
Definition: Style.C:21
Double_t y1[nxs]
Definition: Style.C:20
Double_t x2[nxs]
Definition: Style.C:19
G4double G4NeutronHPJENDLHEData::G4double result
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
Double_t y
Definition: plot.C:279
G4double GetWeightedBinIntegral(const G4InterpolationScheme &aScheme, const G4double x1, const G4double x2, const G4double y1, const G4double y2)
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
G4double GetBinIntegral(const G4InterpolationScheme &aScheme, const G4double x1, const G4double x2, const G4double y1, const G4double y2)
Double_t x1[nxs]
Definition: Style.C:18
G4InterpolationScheme
#define G4endl
Definition: G4ios.hh:61
G4double Lin(G4double x, G4double x1, G4double x2, G4double y1, G4double y2)
double G4double
Definition: G4Types.hh:76