| Geant4 Cross Reference |
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 // neutron_hp -- source file 26 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996 27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans 28 // A prototype of the low energy neutron transport model.
29 // 29 //
30 // P. Arce, June-2014 Conversion neutron_hp to 30 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
31 // 31 //
32 #include "G4ParticleHPInterpolator.hh" 32 #include "G4ParticleHPInterpolator.hh"
33 <<
34 #include "G4Pow.hh" 33 #include "G4Pow.hh"
35 34
36 G4double G4ParticleHPInterpolator::GetBinInteg << 35 G4double G4ParticleHPInterpolator::
37 << 36 GetBinIntegral(const G4InterpolationScheme & aScheme,
38 << 37 const G4double x1,const G4double x2,const G4double y1,const G4double y2)
39 { // inline again later on @@@@ << 38 { // inline again later on @@@@
40 G4double result = 0; << 39 G4double result = 0;
41 if (aScheme == HISTO || aScheme == CHISTO || << 40 if(aScheme==HISTO||aScheme==CHISTO||aScheme==UHISTO)
42 result = y1 * (x2 - x1); << 41 {
43 } << 42 result = y1*(x2-x1);
44 else if (aScheme == LINLIN || aScheme == CLI << 43 }
45 result = 0.5 * (y2 + y1) * (x2 - x1); << 44 else if(aScheme==LINLIN||aScheme==CLINLIN||aScheme==ULINLIN)
46 } << 45 {
47 else if (aScheme == LINLOG || aScheme == CLI << 46 result = 0.5*(y2+y1)*(x2-x1);
48 if (x1 == 0) << 47 }
49 result = y1; << 48 else if(aScheme==LINLOG||aScheme==CLINLOG||aScheme==ULINLOG)
50 else if (x2 == 0) << 49 {
51 result = y2; << 50 if(x1==0) result = y1;
52 else { << 51 else if(x2==0) result = y2;
53 G4double b = (y2 - y1) / (G4Log(x2) - G4 << 52 else
54 G4double a = y1 - b * G4Log(x1); << 53 {
55 result = (a - b) * (x2 - x1) + b * (x2 * << 54 G4double b = (y2-y1)/(G4Log(x2)-G4Log(x1));
>> 55 G4double a = y1 - b*G4Log(x1);
>> 56 result = (a-b)*(x2-x1) + b*(x2*G4Log(x2)-x1*G4Log(x1));
>> 57 }
56 } 58 }
57 } << 59 else if(aScheme==LOGLIN||aScheme==CLOGLIN||aScheme==ULOGLIN)
58 else if (aScheme == LOGLIN || aScheme == CLO << 60 {
59 if (y1 == 0 || y2 == 0) { << 61 if ( y1==0||y2==0 ) {
60 result = 0; << 62 result =0;
61 } << 63 } else {
62 else { << 64 // G4double b = (std::log(y2)-std::log(y1))/(x2-x1);
63 // G4double b = (std::log(y2)-std::log(y << 65 // G4double a = std::log(y1) - b*x1;
64 // G4double a = std::log(y1) - b*x1; << 66 //***************************************************************
65 //************************************** << 67 //EMendoza:
66 // EMendoza: << 68 //result = (std::exp(a)/b)*(std::exp(b*x2)-std::exp(b*x1));
67 // result = (std::exp(a)/b)*(std::exp(b* << 69 //***************************************************************
68 //************************************** << 70 if ( y1!=y2 ) {
69 if (y1 != y2) { << 71 result = (x2-x1)*(y2-y1)/G4Log(y2/y1);
70 result = (x2 - x1) * (y2 - y1) / G4Log << 72 } else {
>> 73 result = y2*(x2-x1);
>> 74 }
>> 75 //***************************************************************
71 } 76 }
72 else { << 77 }
73 result = y2 * (x2 - x1); << 78 else if(aScheme==LOGLOG||aScheme==CLOGLOG||aScheme==ULOGLOG)
>> 79 {
>> 80 if(x1==0) result = y1;
>> 81 else if(x2==0) result = y2;
>> 82 else if(y1==0||y2==0) result =0;
>> 83 else
>> 84 {
>> 85 G4double b = (G4Log(y2)-G4Log(y1))/(G4Log(x2)-G4Log(x1));
>> 86 G4double a = G4Log(y1) - b*G4Log(x1);;
>> 87 result = (G4Exp(a)/(b+1))*(G4Pow::GetInstance()->powA(x2,b+1)-G4Pow::GetInstance()->powA(x1,b+1));
74 } 88 }
75 //************************************** <<
76 } 89 }
77 } << 90 else
78 else if (aScheme == LOGLOG || aScheme == CLO << 91 {
79 if (x1 == 0) << 92 throw G4HadronicException(__FILE__, __LINE__, "Unknown interpolation scheme in G4ParticleHPVector::Integrate");
80 result = y1; << 93 }
81 else if (x2 == 0) << 94 return result;
82 result = y2; << 95 }
83 else if (y1 == 0 || y2 == 0) << 96 G4double G4ParticleHPInterpolator::
84 result = 0; << 97 GetWeightedBinIntegral(const G4InterpolationScheme & aScheme,
85 else { << 98 const G4double x1,const G4double x2,const G4double y1,const G4double y2)
86 G4double b = (G4Log(y2) - G4Log(y1)) / ( << 99 { // inline again later on @@@@
87 G4double a = G4Log(y1) - b * G4Log(x1); << 100 G4double result = 0;
88 ; << 101 if(aScheme==HISTO||aScheme==CHISTO||aScheme==UHISTO)
89 result = (G4Exp(a) / (b + 1)) << 102 {
90 * (G4Pow::GetInstance()->powA(x << 103 result = 0.5*y1*(x2*x2-x1*x1);
>> 104 }
>> 105 else if(aScheme==LINLIN||aScheme==CLINLIN||aScheme==ULINLIN)
>> 106 {
>> 107 // G4double b = (y2-y1)/(x2-x1);
>> 108 // G4double a = y1 - b*x1;
>> 109 // result = 0.5*a*(x2*x2-x1*x1) + (b/3.)*(x2*x2*x2-x1*x1*x1);
>> 110 // Factor out x2-x1 to avoid divide by zero
>> 111
>> 112 result = (y1*x2 - y2*x1)*(x2 + x1)/2. + (y2-y1)*(x2*x2 + x2*x1 + x1*x1)/3.;
>> 113 }
>> 114 else if(aScheme==LINLOG||aScheme==CLINLOG||aScheme==ULINLOG)
>> 115 {
>> 116 if(x1==0) result = y1;
>> 117 else if(x2==0) result = y2;
>> 118 else
>> 119 {
>> 120 G4double b = (y2-y1)/(G4Log(x2)-G4Log(x1));
>> 121 G4double a = y1 - b*G4Log(x1);
>> 122 result = ( x2*x2/2. * (a-b/2.+b*G4Log(x2)) )
>> 123 -( x1*x1/2. * (a-b/2.+b*G4Log(x1)) );
>> 124 }
91 } 125 }
92 } << 126 else if(aScheme==LOGLIN||aScheme==CLOGLIN||aScheme==ULOGLIN)
93 else { << 127 {
94 throw G4HadronicException(__FILE__, __LINE << 128 if(y1==0||y2==0) result = 0;
95 "Unknown interpo << 129 else
96 } << 130 {
97 return result; << 131 G4double b = (G4Log(y2)-G4Log(y1))/(x2-x1);
98 } << 132 G4double a = G4Log(y1) - b*x1;
99 G4double G4ParticleHPInterpolator::GetWeighted << 133 result = G4Exp(a)/(b*b)*( G4Exp(b*x2)*(b*x2-1.) - G4Exp(b*x1)*(b*x1-1.) );
100 << 134 }
101 <<
102 { // inline again later on @@@@ <<
103 G4double result = 0; <<
104 if (aScheme == HISTO || aScheme == CHISTO || <<
105 result = 0.5 * y1 * (x2 * x2 - x1 * x1); <<
106 } <<
107 else if (aScheme == LINLIN || aScheme == CLI <<
108 // G4double b = (y2-y1)/(x2-x1); <<
109 // G4double a = y1 - b*x1; <<
110 // result = 0.5*a*(x2*x2-x1*x1) + ( <<
111 // Factor out x2-x1 to avoid divide by ze <<
112 <<
113 result = (y1 * x2 - y2 * x1) * (x2 + x1) / <<
114 } <<
115 else if (aScheme == LINLOG || aScheme == CLI <<
116 if (x1 == 0) <<
117 result = y1; <<
118 else if (x2 == 0) <<
119 result = y2; <<
120 else { <<
121 G4double b = (y2 - y1) / (G4Log(x2) - G4 <<
122 G4double a = y1 - b * G4Log(x1); <<
123 result = (x2 * x2 / 2. * (a - b / 2. + b <<
124 - (x1 * x1 / 2. * (a - b / 2. + <<
125 } 135 }
126 } << 136 else if(aScheme==LOGLOG||aScheme==CLOGLOG||aScheme==ULOGLOG)
127 else if (aScheme == LOGLIN || aScheme == CLO << 137 {
128 if (y1 == 0 || y2 == 0) << 138 if(x1==0) result = y1;
129 result = 0; << 139 else if(x2==0) result = y2;
130 else { << 140 else if(y1==0||y2==0) result = 0;
131 G4double b = (G4Log(y2) - G4Log(y1)) / ( << 141 else
132 G4double a = G4Log(y1) - b * x1; << 142 {
133 result = G4Exp(a) / (b * b) * (G4Exp(b * << 143 G4double b = (G4Log(y2)-G4Log(y1))/(G4Log(x2)-G4Log(x1));
>> 144 G4double a = G4Log(y1) - b*G4Log(x1);;
>> 145 result = G4Exp(a)/(b+2.)*( G4Pow::GetInstance()->powA(x2, b+2.) - G4Pow::GetInstance()->powA(x1, b+2) );
>> 146 }
134 } 147 }
135 } << 148 else
136 else if (aScheme == LOGLOG || aScheme == CLO << 149 {
137 if (x1 == 0) << 150 throw G4HadronicException(__FILE__, __LINE__, "Unknown interpolation scheme in G4ParticleHPVector::Integrate");
138 result = y1; <<
139 else if (x2 == 0) <<
140 result = y2; <<
141 else if (y1 == 0 || y2 == 0) <<
142 result = 0; <<
143 else { <<
144 G4double b = (G4Log(y2) - G4Log(y1)) / ( <<
145 G4double a = G4Log(y1) - b * G4Log(x1); <<
146 ; <<
147 result = G4Exp(a) / (b + 2.) <<
148 * (G4Pow::GetInstance()->powA(x <<
149 } 151 }
>> 152 return result;
150 } 153 }
151 else { <<
152 throw G4HadronicException(__FILE__, __LINE <<
153 "Unknown interpo <<
154 } <<
155 return result; <<
156 } <<
157 154