Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4ParticleHPLegendreStore Class Reference

#include <G4ParticleHPLegendreStore.hh>

Public Member Functions

 G4ParticleHPLegendreStore (G4int n)
 
 ~G4ParticleHPLegendreStore ()
 
void Init (G4int i, G4double e, G4int n)
 
void SetNPoints (G4int n)
 
void SetEnergy (G4int i, G4double energy)
 
void SetTemperature (G4int i, G4double temp)
 
void SetCoeff (G4int i, G4int l, G4double coeff)
 
void SetCoeff (G4int i, G4ParticleHPLegendreTable *theTable)
 
G4double GetCoeff (G4int i, G4int l)
 
G4double GetEnergy (G4int i)
 
G4double GetTemperature (G4int i)
 
G4int GetNumberOfPoly (G4int i)
 
G4double SampleDiscreteTwoBody (G4double anEnergy)
 
G4double SampleElastic (G4double anEnergy)
 
G4double Sample (G4double energy)
 
G4double SampleMax (G4double energy)
 
G4double Integrate (G4int k, G4double costh)
 
void InitInterpolation (std::istream &aDataFile)
 
void SetManager (G4InterpolationManager &aManager)
 

Detailed Description

Definition at line 39 of file G4ParticleHPLegendreStore.hh.

Constructor & Destructor Documentation

◆ G4ParticleHPLegendreStore()

G4ParticleHPLegendreStore::G4ParticleHPLegendreStore ( G4int n)
inline

Definition at line 42 of file G4ParticleHPLegendreStore.hh.

43 {
44 theCoeff.resize(n);
45 nEnergy = n;
46 }

◆ ~G4ParticleHPLegendreStore()

G4ParticleHPLegendreStore::~G4ParticleHPLegendreStore ( )
inline

Definition at line 48 of file G4ParticleHPLegendreStore.hh.

48{}

Member Function Documentation

◆ GetCoeff()

G4double G4ParticleHPLegendreStore::GetCoeff ( G4int i,
G4int l )
inline

Definition at line 64 of file G4ParticleHPLegendreStore.hh.

64{ return theCoeff[i].GetCoeff(l); }

Referenced by SampleDiscreteTwoBody(), SampleElastic(), and SampleMax().

◆ GetEnergy()

G4double G4ParticleHPLegendreStore::GetEnergy ( G4int i)
inline

Definition at line 65 of file G4ParticleHPLegendreStore.hh.

65{ return theCoeff[i].GetEnergy(); }

Referenced by Sample(), SampleDiscreteTwoBody(), SampleElastic(), and SampleMax().

◆ GetNumberOfPoly()

G4int G4ParticleHPLegendreStore::GetNumberOfPoly ( G4int i)
inline

Definition at line 67 of file G4ParticleHPLegendreStore.hh.

67{ return theCoeff[i].GetNumberOfPoly(); }

◆ GetTemperature()

G4double G4ParticleHPLegendreStore::GetTemperature ( G4int i)
inline

Definition at line 66 of file G4ParticleHPLegendreStore.hh.

66{ return theCoeff[i].GetTemperature(); }

◆ Init()

void G4ParticleHPLegendreStore::Init ( G4int i,
G4double e,
G4int n )
inline

Definition at line 50 of file G4ParticleHPLegendreStore.hh.

50{ theCoeff[i].Init(e, n); }

Referenced by G4ParticleHPAngular::Init(), G4ParticleHPElasticFS::Init(), and G4ParticleHPContAngularPar::Sample().

◆ InitInterpolation()

void G4ParticleHPLegendreStore::InitInterpolation ( std::istream & aDataFile)
inline

Definition at line 75 of file G4ParticleHPLegendreStore.hh.

75{ theManager.Init(aDataFile); }
void Init(G4int aScheme, G4int aRange)

Referenced by G4ParticleHPAngular::Init(), and G4ParticleHPElasticFS::Init().

◆ Integrate()

G4double G4ParticleHPLegendreStore::Integrate ( G4int k,
G4double costh )

Definition at line 301 of file G4ParticleHPLegendreStore.cc.

303{
304 G4double result = 0.;
306 // G4cout <<"the COEFFS "<<k<<" ";
307 // G4cout <<theCoeff[k].GetNumberOfPoly()<<" ";
308 for (G4int l = 0; l < theCoeff[k].GetNumberOfPoly(); l++) {
309 result += theCoeff[k].GetCoeff(l) * theLeg.Integrate(l, costh);
310 // G4cout << theCoeff[k].GetCoeff(l)<<" ";
311 }
312 // G4cout <<G4endl;
313 return result;
314}
double G4double
Definition G4Types.hh:83
int G4int
Definition G4Types.hh:85
G4double Integrate(G4int l, G4double costh)

Referenced by Sample().

◆ Sample()

G4double G4ParticleHPLegendreStore::Sample ( G4double energy)

Definition at line 250 of file G4ParticleHPLegendreStore.cc.

251{
252 G4int i0;
253 G4int low(0), high(0);
254 // G4cout << "G4ParticleHPLegendreStore::Sample "<<energy<<" "<<energy<<" "<<nEnergy<<G4endl;
255 for (i0 = 0; i0 < nEnergy; i0++) {
256 // G4cout <<"theCoeff["<<i0<<"].GetEnergy() = "<<theCoeff[i0].GetEnergy()<<G4endl;
257 high = i0;
258 if (theCoeff[i0].GetEnergy() > energy) break;
259 }
260 low = std::max(0, high - 1);
261 // G4cout << "G4ParticleHPLegendreStore::Sample high, low: "<<high<<", "<<low<<G4endl;
262 G4ParticleHPVector theBuffer;
264 G4double x1, x2, y1, y2, y;
265 x1 = theCoeff[low].GetEnergy();
266 x2 = theCoeff[high].GetEnergy();
267 // G4cout << "the xes "<<x1<<" "<<x2<<G4endl;
268 G4double costh = 0;
269 for (i0 = 0; i0 < 601; i0++) {
270 costh = G4double(i0 - 300) / 300.;
271 y1 = Integrate(low, costh);
272 y2 = Integrate(high, costh);
273 y = theInt.Interpolate(theManager.GetScheme(high), energy, x1, x2, y1, y2);
274 theBuffer.SetData(i0, costh, y);
275 // G4cout << "Integration "<<low<<" "<<costh<<" "<<y1<<" "<<y2<<" "<<y<<G4endl;
276 }
277 G4double rand = G4UniformRand();
278 G4int it;
279 for (i0 = 1; i0 < 601; i0++) {
280 it = i0;
281 if (rand < theBuffer.GetY(i0) / theBuffer.GetY(600)) break;
282 // G4cout <<"sampling now "<<i0<<" "
283 // << theBuffer.GetY(i0)<<" "
284 // << theBuffer.GetY(600)<<" "
285 // << rand<<" "
286 // << theBuffer.GetY(i0)/theBuffer.GetY(600)<<G4endl;;
287 }
288 if (it == 601) it = 600;
289 // G4cout << "G4ParticleHPLegendreStore::Sample it "<<rand<<" "<<it<<G4endl;
290 G4double norm = theBuffer.GetY(600);
291 if (norm == 0) return -DBL_MAX;
292 x1 = theBuffer.GetY(it) / norm;
293 x2 = theBuffer.GetY(it - 1) / norm;
294 y1 = theBuffer.GetX(it);
295 y2 = theBuffer.GetX(it - 1);
296 // G4cout << "G4ParticleHPLegendreStore::Sample x y "<<x1<<" "<<y1<<" "<<x2<<" "<<y2<<G4endl;
297 return theInt.Interpolate(theManager.GetScheme(high), rand, x1, x2, y1, y2);
298}
#define G4UniformRand()
Definition Randomize.hh:52
G4InterpolationScheme GetScheme(G4int index) const
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
G4double Integrate(G4int k, G4double costh)
void SetData(G4int i, G4double x, G4double y)
G4double GetY(G4double x)
G4double GetX(G4int i) const
#define DBL_MAX
Definition templates.hh:62

◆ SampleDiscreteTwoBody()

G4double G4ParticleHPLegendreStore::SampleDiscreteTwoBody ( G4double anEnergy)

Definition at line 45 of file G4ParticleHPLegendreStore.cc.

46{
47 G4double result = 0.;
48
49 G4int i0;
50 G4int low(0), high(0);
52 for (i0 = 0; i0 < nEnergy; i0++) {
53 high = i0;
54 if (theCoeff[i0].GetEnergy() > anEnergy) break;
55 }
56 low = std::max(0, high - 1);
58 G4double x, x1, x2;
59 x = anEnergy;
60 x1 = theCoeff[low].GetEnergy();
61 x2 = theCoeff[high].GetEnergy();
62 G4double theNorm = 0;
63 G4double try01 = 0, try02 = 0;
64 G4double max1, max2, costh;
65 max1 = 0;
66 max2 = 0;
67 G4int l, m_tmp;
68 for (i0 = 0; i0 < 601; i0++) {
69 costh = G4double(i0 - 300) / 300.;
70 try01 = 0.5;
71 for (m_tmp = 0; m_tmp < theCoeff[low].GetNumberOfPoly(); m_tmp++) {
72 l = m_tmp + 1;
73 try01 += (2. * l + 1) / 2. * theCoeff[low].GetCoeff(m_tmp) * theLeg.Evaluate(l, costh);
74 }
75 if (try01 > max1) max1 = try01;
76 try02 = 0.5;
77 for (m_tmp = 0; m_tmp < theCoeff[high].GetNumberOfPoly(); m_tmp++) {
78 l = m_tmp + 1;
79 try02 += (2. * l + 1) / 2. * theCoeff[high].GetCoeff(m_tmp) * theLeg.Evaluate(l, costh);
80 }
81 if (try02 > max2) max2 = try02;
82 }
83 theNorm = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, max1, max2);
84
85 G4double value, random;
86 G4double v1, v2;
87 G4int icounter = 0;
88 G4int icounter_max = 1024;
89 do {
90 icounter++;
91 if (icounter > icounter_max) {
92 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of "
93 << __FILE__ << "." << G4endl;
94 break;
95 }
96 v1 = 0.5;
97 v2 = 0.5;
98 result = 2. * G4UniformRand() - 1.;
99 for (m_tmp = 0; m_tmp < theCoeff[low].GetNumberOfPoly(); m_tmp++) {
100 l = m_tmp + 1;
101 G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN
102 v1 += (2. * l + 1) / 2. * theCoeff[low].GetCoeff(m_tmp) * legend;
103 }
104 for (m_tmp = 0; m_tmp < theCoeff[high].GetNumberOfPoly(); m_tmp++) {
105 l = m_tmp + 1;
106 G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN
107 v2 += (2. * l + 1) / 2. * theCoeff[high].GetCoeff(m_tmp) * legend;
108 }
109 // v1 = std::max(0.,v1); // Workaround in case one of the distributions is fully non-physical.
110 // v2 = std::max(0.,v2);
111 value = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, v1, v2);
112 random = G4UniformRand();
113 if (0 >= theNorm) break; // Workaround for negative cross-section values. @@@@ 31 May 2000
114 } while (random > value / theNorm); // Loop checking, 11.05.2015, T. Koi
115
116 return result;
117}
#define G4endl
Definition G4ios.hh:67
G4GLOB_DLL std::ostream G4cout
G4double Evaluate(G4int l, G4double costh)
G4double GetCoeff(G4int i, G4int l)

Referenced by G4ParticleHPDiscreteTwoBody::Sample().

◆ SampleElastic()

G4double G4ParticleHPLegendreStore::SampleElastic ( G4double anEnergy)

Definition at line 188 of file G4ParticleHPLegendreStore.cc.

189{
190 G4double result = 0.;
191
192 G4int i0;
193 G4int low(0), high(0);
195 for (i0 = 0; i0 < nEnergy; i0++) {
196 high = i0;
197 if (theCoeff[i0].GetEnergy() > anEnergy) break;
198 }
199 low = std::max(0, high - 1);
201 G4double x, x1, x2;
202 x = anEnergy;
203 x1 = theCoeff[low].GetEnergy();
204 x2 = theCoeff[high].GetEnergy();
205 G4double theNorm = 0;
206 G4double try01 = 0, try02 = 0, try11 = 0, try12 = 0;
207 G4double try1, try2;
208 G4int l;
209 for (l = 0; l < theCoeff[low].GetNumberOfPoly(); l++) {
210 try01 += (2. * l + 1) / 2. * theCoeff[low].GetCoeff(l) * theLeg.Evaluate(l, -1.);
211 try11 += (2. * l + 1) / 2. * theCoeff[low].GetCoeff(l) * theLeg.Evaluate(l, +1.);
212 }
213 for (l = 0; l < theCoeff[high].GetNumberOfPoly(); l++) {
214 try02 += (2. * l + 1) / 2. * theCoeff[high].GetCoeff(l) * theLeg.Evaluate(l, -1.);
215 try12 += (2. * l + 1) / 2. * theCoeff[high].GetCoeff(l) * theLeg.Evaluate(l, +1.);
216 }
217 try1 = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, try01, try02);
218 try2 = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, try11, try12);
219 theNorm = std::max(try1, try2);
220
221 G4double value, random;
222 G4double v1, v2;
223 G4int icounter = 0;
224 G4int icounter_max = 1024;
225 do {
226 icounter++;
227 if (icounter > icounter_max) {
228 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of "
229 << __FILE__ << "." << G4endl;
230 break;
231 }
232 v1 = 0;
233 v2 = 0;
234 result = 2. * G4UniformRand() - 1.;
235 for (l = 0; l < theCoeff[low].GetNumberOfPoly(); l++) {
236 G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN
237 v1 += (2. * l + 1) / 2. * theCoeff[low].GetCoeff(l) * legend;
238 }
239 for (l = 0; l < theCoeff[high].GetNumberOfPoly(); l++) {
240 G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN
241 v2 += (2. * l + 1) / 2. * theCoeff[high].GetCoeff(l) * legend;
242 }
243 value = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, v1, v2);
244 random = G4UniformRand();
245 } while (random > value / theNorm); // Loop checking, 11.05.2015, T. Koi
246
247 return result;
248}

Referenced by G4ParticleHPElasticFS::ApplyYourself().

◆ SampleMax()

G4double G4ParticleHPLegendreStore::SampleMax ( G4double energy)

Definition at line 119 of file G4ParticleHPLegendreStore.cc.

120{
121 G4double result = 0.;
122
123 G4int i0;
124 G4int low(0), high(0);
126 for (i0 = 0; i0 < nEnergy; i0++) {
127 high = i0;
128 if (theCoeff[i0].GetEnergy() > anEnergy) break;
129 }
130 low = std::max(0, high - 1);
132 G4double x, x1, x2;
133 x = anEnergy;
134 x1 = theCoeff[low].GetEnergy();
135 x2 = theCoeff[high].GetEnergy();
136 G4double theNorm = 0;
137 G4double try01 = 0, try02 = 0;
138 G4double max1, max2, costh;
139 max1 = 0;
140 max2 = 0;
141 G4int l;
142 for (i0 = 0; i0 < 601; i0++) {
143 costh = G4double(i0 - 300) / 300.;
144 try01 = 0;
145 for (l = 0; l < theCoeff[low].GetNumberOfPoly(); l++) {
146 try01 += (2. * l + 1) / 2. * theCoeff[low].GetCoeff(l) * theLeg.Evaluate(l, costh);
147 }
148 if (try01 > max1) max1 = try01;
149 try02 = 0;
150 for (l = 0; l < theCoeff[high].GetNumberOfPoly(); l++) {
151 try02 += (2. * l + 1) / 2. * theCoeff[high].GetCoeff(l) * theLeg.Evaluate(l, costh);
152 }
153 if (try02 > max2) max2 = try02;
154 }
155 theNorm = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, max1, max2);
156
157 G4double value, random;
158 G4double v1, v2;
159 G4int icounter = 0;
160 G4int icounter_max = 1024;
161 do {
162 icounter++;
163 if (icounter > icounter_max) {
164 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of "
165 << __FILE__ << "." << G4endl;
166 break;
167 }
168 v1 = 0;
169 v2 = 0;
170 result = 2. * G4UniformRand() - 1.;
171 for (l = 0; l < theCoeff[low].GetNumberOfPoly(); l++) {
172 G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN
173 v1 += (2. * l + 1) / 2. * theCoeff[low].GetCoeff(l) * legend;
174 }
175 for (l = 0; l < theCoeff[high].GetNumberOfPoly(); l++) {
176 G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN
177 v2 += (2. * l + 1) / 2. * theCoeff[high].GetCoeff(l) * legend;
178 }
179 v1 = std::max(0., v1); // Workaround in case one of the distributions is fully non-physical.
180 v2 = std::max(0., v2);
181 value = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, v1, v2);
182 random = G4UniformRand();
183 if (0 >= theNorm) break; // Workaround for negative cross-section values. @@@@ 31 May 2000
184 } while (random > value / theNorm); // Loop checking, 11.05.2015, T. Koi
185 return result;
186}

Referenced by G4ParticleHPPhotonDist::GetPhotons(), G4ParticleHPContAngularPar::Sample(), and G4ParticleHPAngular::SampleAndUpdate().

◆ SetCoeff() [1/2]

void G4ParticleHPLegendreStore::SetCoeff ( G4int i,
G4int l,
G4double coeff )
inline

◆ SetCoeff() [2/2]

void G4ParticleHPLegendreStore::SetCoeff ( G4int i,
G4ParticleHPLegendreTable * theTable )
inline

Definition at line 55 of file G4ParticleHPLegendreStore.hh.

56 {
57 if (i > nEnergy)
58 throw G4HadronicException(__FILE__, __LINE__, "LegendreTableIndex out of range");
59 theCoeff[i] = *theTable;
60 // not here -- see G4ParticleHPPhotonDist.cc line 275
61 // delete theTable;
62 }

◆ SetEnergy()

void G4ParticleHPLegendreStore::SetEnergy ( G4int i,
G4double energy )
inline

Definition at line 52 of file G4ParticleHPLegendreStore.hh.

52{ theCoeff[i].SetEnergy(energy); }

◆ SetManager()

void G4ParticleHPLegendreStore::SetManager ( G4InterpolationManager & aManager)
inline

Definition at line 77 of file G4ParticleHPLegendreStore.hh.

77{ theManager = aManager; }

Referenced by G4ParticleHPContAngularPar::Sample(), and G4ParticleHPDiscreteTwoBody::Sample().

◆ SetNPoints()

void G4ParticleHPLegendreStore::SetNPoints ( G4int n)
inline

Definition at line 51 of file G4ParticleHPLegendreStore.hh.

51{ nEnergy = n; }

◆ SetTemperature()

void G4ParticleHPLegendreStore::SetTemperature ( G4int i,
G4double temp )
inline

Definition at line 53 of file G4ParticleHPLegendreStore.hh.

53{ theCoeff[i].SetTemperature(temp); }

Referenced by G4ParticleHPAngular::Init(), and G4ParticleHPElasticFS::Init().


The documentation for this class was generated from the following files: