ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNACPA100LogLogInterpolation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNACPA100LogLogInterpolation.cc
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 // Based on the work of M. Terrissol and M. C. Bordage
27 //
28 // Users are requested to cite the following papers:
29 // - M. Terrissol, A. Baudre, Radiat. Prot. Dosim. 31 (1990) 175-177
30 // - M.C. Bordage, J. Bordes, S. Edel, M. Terrissol, X. Franceries,
31 // M. Bardies, N. Lampe, S. Incerti, Phys. Med. 32 (2016) 1833-1840
32 //
33 // Authors of this class:
34 // M.C. Bordage, M. Terrissol, S. Edel, J. Bordes, S. Incerti
35 //
36 // 15.01.2014: creation
37 //
38 
40 
41 // Constructor
42 
44 { }
45 
46 // Destructor
47 
49 { }
50 
52 { return new G4DNACPA100LogLogInterpolation; }
53 
54 
56  const G4DataVector& points,
57  const G4DataVector& data) const
58 {
59  //G4cout << "G4DNACPA100LogLogInterpolation is performed (2 arguments) " << G4endl;
60  G4int nBins = data.size() - 1;
61 //G4double oldresult = 0.;
62  G4double value = 0.;
63  if (x < points[0])
64  {
65  value = 0.;
66  }
67  else if (bin < nBins)
68  {
69  G4double e1 = points[bin];
70  G4double e2 = points[bin+1];
71  G4double d1 = data[bin];
72  G4double d2 = data[bin+1];
73 // Check of e1, e2, d1 and d2 values to avoid floating-point errors when estimating the interpolated value below -- S.I., Jun. 2008
74  if ((d1 > 0.) && (d2 > 0.) && (e1 > 0.) && (e2 > 0.))
75  {
76 // Streamline the Log-Log Interpolation formula in order to reduce the required number of log10() function calls
77 // Variable oldresult contains the result of old implementation of Log-Log interpolation -- M.G.P. Jun. 2001
78 // oldresult = (std::log10(d1)*std::log10(e2/x) + std::log10(d2)*std::log10(x/e1)) / std::log10(e2/e1);
79 // oldresult = std::pow(10.,oldresult);
80 // Variable value contains the result of new implementation, after streamlining the math operation -- N.A.K. Oct. 2008
81  value = std::log10(d1)+(std::log10(d2/d1)/std::log10(e2/e1)*std::log10(x/e1));
82  value = std::pow(10.,value);
83 // Test of the new implementation result (value variable) against the old one (oldresult) -- N.A.K. Dec. 2008
84 // G4double diffResult = value - oldresult;
85 // G4double relativeDiff = 1e-11;
86 // Comparison of the two values based on a max allowable relative difference
87 // if ( std::fabs(diffResult) > relativeDiff*std::fabs(oldresult) )
88 // {
89 // Abort comparison when at least one of two results is infinite
90 // if ((!std::isinf(oldresult)) && (!std::isinf(value)))
91 // {
92 // G4cout << "G4DNACPA100LogLogInterpolation> Old Interpolated Value is:" << oldresult << G4endl;
93 // G4cout << "G4DNACPA100LogLogInterpolation> New Interpolated Value is:" << value << G4endl << G4endl;
94 // G4cerr << "G4DNACPA100LogLogInterpolation> Error in Interpolation:" << G4endl;
95 // G4cerr << "The difference between new and old interpolated value is:" << diffResult << G4endl << G4endl;
96 // }
97 // }
98  }
99  else value = 0.;
100  }
101  else
102  {
103  value = data[nBins];
104  }
105  return value;
106 }
107 
108 
109 // Nicolas A. Karakatsanis: New implementation of log-log interpolation after directly loading
110 // logarithmic values from G4EMLOW dataset
111 
113  const G4DataVector& points,
114  const G4DataVector& data,
115  const G4DataVector& log_points,
116  const G4DataVector& log_data) const
117 {
118  G4int nBins = data.size() - 1;
119  G4double value = 0.;
120  G4double log_x = std::log10(x);
121  if (x < points[0])
122  {
123  value = 0.;
124  }
125  else if (bin < nBins)
126  {
127  G4double log_e1 = log_points[bin];
128  //G4double log_e2 = log_points[bin+1];
129  G4double log_d1 = log_data[bin];
130  G4double log_d2 = log_data[bin+1];
131 
132  //G4cout << "x = " << x << " , logx = " << log_x << " , bin = " << bin << G4endl;
133  //G4cout << "e1 = " << points[bin] << " d1 = " << data[bin] << G4endl;
134  //G4cout << "e2 = " << points[bin+1] << " d2 = " << data[bin+1] << G4endl;
135  //G4cout << "loge1 = " << log_e1 << " logd1 = " << log_d1 << G4endl;
136  //G4cout << "loge2 = " << log_e2 << " logd2 = " << log_d2 << G4endl;
137  //G4cout << "interpol " << log_d1 + (log_d2 - log_d1)*(log_x - log_e1)/(log_e2 - log_e1) << " " << G4endl;
138 
139 
140  // Values e1, e2, d1 and d2 are the log values of the corresponding
141  // original energy and data values. Simple linear interpolation performed
142  // on loagarithmic data should be equivalent to log-log interpolation
143 
144  // CPA100 specific
145  //value = log_d1 + (log_d2 - log_d1)*(log_x - log_e1)/(log_e2 - log_e1);
146  // value = log_d1;
147 
148  value = log_d2; // UPPER VALUE INTERPOLATION
149  if (log_x == log_e1) value = log_d1; // IN CASE OF EQUALITY
150 
151 // Delogarithmize to obtain interpolated value
152  value = std::pow(10.,value);
153  }
154  else
155  {
156  value = data[nBins];
157  }
158 
159  return value;
160 }