Geant4_10
G4NeutronHPInterpolator.cc
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 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron transport model.
29 //
31 
34  const G4double x1,const G4double x2,const G4double y1,const G4double y2)
35  { // inline again later on @@@@
36  G4double result = 0;
37  if(aScheme==HISTO||aScheme==CHISTO||aScheme==UHISTO)
38  {
39  result = y1*(x2-x1);
40  }
41  else if(aScheme==LINLIN||aScheme==CLINLIN||aScheme==ULINLIN)
42  {
43  result = 0.5*(y2+y1)*(x2-x1);
44  }
45  else if(aScheme==LINLOG||aScheme==CLINLOG||aScheme==ULINLOG)
46  {
47  if(x1==0) result = y1;
48  else if(x2==0) result = y2;
49  else
50  {
51  G4double b = (y2-y1)/(std::log(x2)-std::log(x1));
52  G4double a = y1 - b*std::log(x1);
53  result = (a-b)*(x2-x1) + b*(x2*std::log(x2)-x1*std::log(x1));
54  }
55  }
56  else if(aScheme==LOGLIN||aScheme==CLOGLIN||aScheme==ULOGLIN)
57  {
58  if(y1==0||y2==0) result =0;
59  else
60  {
61  G4double b = (std::log(y2)-std::log(y1))/(x2-x1);
62  G4double a = std::log(y1) - b*x1;
63  //***************************************************************
64  //EMendoza:
65  //result = (std::exp(a)/b)*(std::exp(b*x2)-std::exp(b*x1));
66  //***************************************************************
67  if(b!=0){
68  result = (std::exp(a)/b)*(std::exp(b*x2)-std::exp(b*x1));
69  }
70  else{
71  result=y2*(x2-x1);
72  }
73  //***************************************************************
74  }
75  }
76  else if(aScheme==LOGLOG||aScheme==CLOGLOG||aScheme==ULOGLOG)
77  {
78  if(x1==0) result = y1;
79  else if(x2==0) result = y2;
80  else if(y1==0||y2==0) result =0;
81  else
82  {
83  G4double b = (std::log(y2)-std::log(y1))/(std::log(x2)-std::log(x1));
84  G4double a = std::log(y1) - b*std::log(x1);;
85  result = (std::exp(a)/(b+1))*(std::pow(x2,b+1)-std::pow(x1,b+1));
86  }
87  }
88  else
89  {
90  throw G4HadronicException(__FILE__, __LINE__, "Unknown interpolation scheme in G4NeutronHPVector::Integrate");
91  }
92  return result;
93  }
96  const G4double x1,const G4double x2,const G4double y1,const G4double y2)
97  { // inline again later on @@@@
98  G4double result = 0;
99  if(aScheme==HISTO||aScheme==CHISTO||aScheme==UHISTO)
100  {
101  result = 0.5*y1*(x2*x2-x1*x1);
102  }
103  else if(aScheme==LINLIN||aScheme==CLINLIN||aScheme==ULINLIN)
104  {
105  // G4double b = (y2-y1)/(x2-x1);
106  // G4double a = y1 - b*x1;
107  // result = 0.5*a*(x2*x2-x1*x1) + (b/3.)*(x2*x2*x2-x1*x1*x1);
108  // Factor out x2-x1 to avoid divide by zero
109 
110  result = (y1*x2 - y2*x1)*(x2 + x1)/2. + (y2-y1)*(x2*x2 + x2*x1 + x1*x1)/3.;
111  }
112  else if(aScheme==LINLOG||aScheme==CLINLOG||aScheme==ULINLOG)
113  {
114  if(x1==0) result = y1;
115  else if(x2==0) result = y2;
116  else
117  {
118  G4double b = (y2-y1)/(std::log(x2)-std::log(x1));
119  G4double a = y1 - b*std::log(x1);
120  result = ( x2*x2/2. * (a-b/2.+b*std::log(x2)) )
121  -( x1*x1/2. * (a-b/2.+b*std::log(x1)) );
122  }
123  }
124  else if(aScheme==LOGLIN||aScheme==CLOGLIN||aScheme==ULOGLIN)
125  {
126  if(y1==0||y2==0) result = 0;
127  else
128  {
129  G4double b = (std::log(y2)-std::log(y1))/(x2-x1);
130  G4double a = std::log(y1) - b*x1;
131  result = std::exp(a)/(b*b)*( std::exp(b*x2)*(b*x2-1.) - std::exp(b*x1)*(b*x1-1.) );
132  }
133  }
134  else if(aScheme==LOGLOG||aScheme==CLOGLOG||aScheme==ULOGLOG)
135  {
136  if(x1==0) result = y1;
137  else if(x2==0) result = y2;
138  if(y1==0||y2==0) result = 0;
139  else
140  {
141  G4double b = (std::log(y2)-std::log(y1))/(std::log(x2)-std::log(x1));
142  G4double a = std::log(y1) - b*std::log(x1);;
143  result = std::exp(a)/(b+2.)*( std::pow(x2, b+2.) - std::pow(x1, b+2) );
144  }
145  }
146  else
147  {
148  throw G4HadronicException(__FILE__, __LINE__, "Unknown interpolation scheme in G4NeutronHPVector::Integrate");
149  }
150  return result;
151  }
Double_t y2[nxs]
Definition: Style.C:21
Double_t y1[nxs]
Definition: Style.C:20
Double_t x2[nxs]
Definition: Style.C:19
tuple a
Definition: test.py:11
G4double G4NeutronHPJENDLHEData::G4double result
tuple b
Definition: test.py:12
G4double GetWeightedBinIntegral(const G4InterpolationScheme &aScheme, const G4double x1, const G4double x2, const G4double y1, const G4double y2)
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
double G4double
Definition: G4Types.hh:76