| Geant4 Cross Reference |
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 // 070606 fix with Valgrind by T. Koi
28 // 080409 Fix div0 error with G4FPE by T. Koi
29 // 080811 Comment out unused method SetBlocked and SetBuffered
30 // Add required cleaning up in CleanUp by T. Koi
31 //
32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
33 //
34 #ifndef G4ParticleHPVector_h
35 #define G4ParticleHPVector_h 1
36
37 #include "G4Exp.hh"
38 #include "G4InterpolationManager.hh"
39 #include "G4Log.hh"
40 #include "G4ParticleHPDataPoint.hh"
41 #include "G4ParticleHPHash.hh"
42 #include "G4ParticleHPInterpolator.hh"
43 #include "G4PhysicsVector.hh"
44 #include "G4Pow.hh"
45 #include "G4ios.hh"
46 #include "Randomize.hh"
47
48 #include <cmath>
49 #include <fstream>
50 #include <vector>
51
52 #if defined WIN32 - VC
53 # include <float.h>
54 #endif
55
56 class G4ParticleHPVector
57 {
58 friend G4ParticleHPVector& operator+(G4ParticleHPVector& left, G4ParticleHPVector& right);
59
60 public:
61 G4ParticleHPVector();
62
63 G4ParticleHPVector(G4int n);
64
65 ~G4ParticleHPVector();
66
67 G4ParticleHPVector& operator=(const G4ParticleHPVector& right);
68
69 inline void SetVerbose(G4int ff) { Verbose = ff; }
70
71 inline void Times(G4double factor)
72 {
73 G4int i;
74 for (i = 0; i < nEntries; i++) {
75 theData[i].SetY(theData[i].GetY() * factor);
76 }
77 if (theIntegral != nullptr) {
78 theIntegral[i] *= factor;
79 }
80 }
81
82 inline void SetPoint(G4int i, const G4ParticleHPDataPoint& it)
83 {
84 G4double x = it.GetX();
85 G4double y = it.GetY();
86 SetData(i, x, y);
87 }
88
89 inline void SetData(G4int i, G4double x, G4double y)
90 {
91 // G4cout <<"G4ParticleHPVector::SetData called"<<nPoints<<" "<<nEntries<<G4endl;
92 Check(i);
93 if (y > maxValue) maxValue = y;
94 theData[i].SetData(x, y);
95 }
96 inline void SetX(G4int i, G4double e)
97 {
98 Check(i);
99 theData[i].SetX(e);
100 }
101 inline void SetEnergy(G4int i, G4double e)
102 {
103 Check(i);
104 theData[i].SetX(e);
105 }
106 inline void SetY(G4int i, G4double x)
107 {
108 Check(i);
109 if (x > maxValue) maxValue = x;
110 theData[i].SetY(x);
111 }
112 inline void SetXsec(G4int i, G4double x)
113 {
114 Check(i);
115 if (x > maxValue) maxValue = x;
116 theData[i].SetY(x);
117 }
118 inline G4double GetEnergy(G4int i) const { return theData[i].GetX(); }
119 inline G4double GetXsec(G4int i) { return theData[i].GetY(); }
120 inline G4double GetX(G4int i) const
121 {
122 if (i < 0) i = 0;
123 if (i >= GetVectorLength()) i = GetVectorLength() - 1;
124 return theData[i].GetX();
125 }
126 inline const G4ParticleHPDataPoint& GetPoint(G4int i) const { return theData[i]; }
127
128 void Hash()
129 {
130 G4int i;
131 G4double x, y;
132 for (i = 0; i < nEntries; i++) {
133 if (0 == (i + 1) % 10) {
134 x = GetX(i);
135 y = GetY(i);
136 theHash.SetData(i, x, y);
137 }
138 }
139 }
140
141 void ReHash()
142 {
143 theHash.Clear();
144 Hash();
145 }
146
147 G4double GetXsec(G4double e);
148
149 G4int GetEnergyIndex(G4double &e); // method added by M. Zmeskal 02/2024
150
151 G4double GetXsec(G4double e, G4int min)
152 {
153 G4int i;
154 for (i = min; i < nEntries; i++) {
155 if (theData[i].GetX() > e) break;
156 }
157 G4int low = i - 1;
158 G4int high = i;
159 if (i == 0) {
160 low = 0;
161 high = 1;
162 }
163 else if (i == nEntries) {
164 low = nEntries - 2;
165 high = nEntries - 1;
166 }
167 G4double y;
168 if (e < theData[nEntries - 1].GetX()) {
169 // Protect against doubled-up x values
170 if ((theData[high].GetX() - theData[low].GetX()) / theData[high].GetX() < 0.000001) {
171 y = theData[low].GetY();
172 }
173 else {
174 y = theInt.Interpolate(theManager.GetScheme(high), e, theData[low].GetX(),
175 theData[high].GetX(), theData[low].GetY(), theData[high].GetY());
176 }
177 }
178 else {
179 y = theData[nEntries - 1].GetY();
180 }
181 return y;
182 }
183
184 inline G4double GetY(G4double x) { return GetXsec(x); }
185 inline G4int GetVectorLength() const { return nEntries; }
186
187 inline G4double GetY(G4int i)
188 {
189 if (i < 0) i = 0;
190 if (i >= GetVectorLength()) i = GetVectorLength() - 1;
191 return theData[i].GetY();
192 }
193
194 inline G4double GetY(G4int i) const
195 {
196 if (i < 0) i = 0;
197 if (i >= GetVectorLength()) i = GetVectorLength() - 1;
198 return theData[i].GetY();
199 }
200 void Dump();
201
202 inline void InitInterpolation(std::istream& aDataFile) { theManager.Init(aDataFile); }
203
204 void Init(std::istream& aDataFile, G4int total, G4double ux = 1., G4double uy = 1.)
205 {
206 G4double x, y;
207 for (G4int i = 0; i < total; i++) {
208 aDataFile >> x >> y;
209 x *= ux;
210 y *= uy;
211 SetData(i, x, y);
212 if (0 == nEntries % 10) {
213 theHash.SetData(nEntries - 1, x, y);
214 }
215 }
216 }
217
218 void Init(std::istream& aDataFile, G4double ux = 1., G4double uy = 1.)
219 {
220 G4int total;
221 aDataFile >> total;
222 delete[] theData;
223 theData = new G4ParticleHPDataPoint[total];
224 nPoints = total;
225 nEntries = 0;
226 theManager.Init(aDataFile);
227 Init(aDataFile, total, ux, uy);
228 }
229
230 void ThinOut(G4double precision);
231
232 inline void SetLabel(G4double aLabel) { label = aLabel; }
233
234 inline G4double GetLabel() { return label; }
235
236 inline void CleanUp()
237 {
238 nEntries = 0;
239 theManager.CleanUp();
240 maxValue = -DBL_MAX;
241 theHash.Clear();
242 // 080811 TK DB
243 delete[] theIntegral;
244 theIntegral = nullptr;
245 }
246
247 // merges the vectors active and passive into *this
248 inline void Merge(G4ParticleHPVector* active, G4ParticleHPVector* passive)
249 {
250 CleanUp();
251 G4int s_tmp = 0, n = 0, m_tmp = 0;
252 G4ParticleHPVector* tmp;
253 G4int a = s_tmp, p = n, t;
254 while (a < active->GetVectorLength()
255 && p < passive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
256 {
257 if (active->GetEnergy(a) <= passive->GetEnergy(p)) {
258 G4double xa = active->GetEnergy(a);
259 G4double yy = active->GetXsec(a);
260 SetData(m_tmp, xa, yy);
261 theManager.AppendScheme(m_tmp, active->GetScheme(a));
262 m_tmp++;
263 a++;
264 G4double xp = passive->GetEnergy(p);
265
266 // 080409 TKDB
267 // if( std::abs(std::abs(xp-xa)/xa)<0.001 ) p++;
268 if (!(xa == 0) && std::abs(std::abs(xp - xa) / xa) < 0.001) p++;
269 }
270 else {
271 tmp = active;
272 t = a;
273 active = passive;
274 a = p;
275 passive = tmp;
276 p = t;
277 }
278 }
279 while (a != active->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
280 {
281 SetData(m_tmp, active->GetEnergy(a), active->GetXsec(a));
282 theManager.AppendScheme(m_tmp++, active->GetScheme(a));
283 a++;
284 }
285 while (p != passive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
286 {
287 if (std::abs(GetEnergy(m_tmp - 1) - passive->GetEnergy(p)) / passive->GetEnergy(p) > 0.001)
288 // if(std::abs(GetEnergy(m)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.001)
289 {
290 SetData(m_tmp, passive->GetEnergy(p), passive->GetXsec(p));
291 theManager.AppendScheme(m_tmp++, active->GetScheme(p));
292 }
293 p++;
294 }
295 }
296
297 void Merge(G4InterpolationScheme aScheme, G4double aValue, G4ParticleHPVector* active,
298 G4ParticleHPVector* passive);
299
300 G4double SampleLin() // Samples X according to distribution Y, linear int
301 {
302 G4double result;
303 if (theIntegral == nullptr) IntegrateAndNormalise();
304 if (GetVectorLength() == 1) {
305 result = theData[0].GetX();
306 }
307 else {
308 G4int i;
309 G4double rand = G4UniformRand();
310
311 // this was replaced
312 // for(i=1;i<GetVectorLength();i++)
313 // {
314 // if(rand<theIntegral[i]/theIntegral[GetVectorLength()-1]) break;
315 // }
316
317 // by this (begin)
318 for (i = GetVectorLength() - 1; i >= 0; i--) {
319 if (rand > theIntegral[i] / theIntegral[GetVectorLength() - 1]) break;
320 }
321 if (i != GetVectorLength() - 1) i++;
322 // until this (end)
323
324 G4double x1, x2, y1, y2;
325 y1 = theData[i - 1].GetX();
326 x1 = theIntegral[i - 1];
327 y2 = theData[i].GetX();
328 x2 = theIntegral[i];
329 if (std::abs((y2 - y1) / y2)
330 < 0.0000001) // not really necessary, since the case is excluded by construction
331 {
332 y1 = theData[i - 2].GetX();
333 x1 = theIntegral[i - 2];
334 }
335 result = theLin.Lin(rand, x1, x2, y1, y2);
336 }
337 return result;
338 }
339
340 G4double Sample(); // Samples X according to distribution Y
341
342 G4double* Debug() { return theIntegral; }
343
344 inline void IntegrateAndNormalise()
345 {
346 G4int i;
347 if (theIntegral != nullptr) return;
348 theIntegral = new G4double[nEntries];
349 if (nEntries == 1) {
350 theIntegral[0] = 1;
351 return;
352 }
353 theIntegral[0] = 0;
354 G4double sum = 0;
355 G4double x1 = 0;
356 G4double x0 = 0;
357 for (i = 1; i < GetVectorLength(); i++) {
358 x1 = theData[i].GetX();
359 x0 = theData[i - 1].GetX();
360 if (std::abs(x1 - x0) > std::abs(x1 * 0.0000001)) {
361 //********************************************************************
362 // EMendoza -> the interpolation scheme is not always lin-lin
363 /*
364 sum+= 0.5*(theData[i].GetY()+theData[i-1].GetY())*
365 (x1-x0);
366 */
367 //********************************************************************
368 G4InterpolationScheme aScheme = theManager.GetScheme(i);
369 G4double y0 = theData[i - 1].GetY();
370 G4double y1 = theData[i].GetY();
371 G4double integ = theInt.GetBinIntegral(aScheme, x0, x1, y0, y1);
372 #if defined WIN32 - VC
373 if (!_finite(integ)) {
374 integ = 0;
375 }
376 #elif defined __IBMCPP__
377 if (isinf(integ) || isnan(integ)) {
378 integ = 0;
379 }
380 #else
381 if (std::isinf(integ) || std::isnan(integ)) {
382 integ = 0;
383 }
384 #endif
385 sum += integ;
386 //********************************************************************
387 }
388 theIntegral[i] = sum;
389 }
390 G4double total = theIntegral[GetVectorLength() - 1];
391 for (i = 1; i < GetVectorLength(); i++) {
392 theIntegral[i] /= total;
393 }
394 }
395
396 inline void Integrate()
397 {
398 G4int i;
399 if (nEntries == 1) {
400 totalIntegral = 0;
401 return;
402 }
403 G4double sum = 0;
404 for (i = 1; i < GetVectorLength(); i++) {
405 if (std::abs((theData[i].GetX() - theData[i - 1].GetX()) / theData[i].GetX()) > 0.0000001) {
406 G4double x1 = theData[i - 1].GetX();
407 G4double x2 = theData[i].GetX();
408 G4double y1 = theData[i - 1].GetY();
409 G4double y2 = theData[i].GetY();
410 G4InterpolationScheme aScheme = theManager.GetScheme(i);
411 if (aScheme == LINLIN || aScheme == CLINLIN || aScheme == ULINLIN) {
412 sum += 0.5 * (y2 + y1) * (x2 - x1);
413 }
414 else if (aScheme == LINLOG || aScheme == CLINLOG || aScheme == ULINLOG) {
415 G4double a = y1;
416 G4double b = (y2 - y1) / (G4Log(x2) - G4Log(x1));
417 sum += (a - b) * (x2 - x1) + b * (x2 * G4Log(x2) - x1 * G4Log(x1));
418 }
419 else if (aScheme == LOGLIN || aScheme == CLOGLIN || aScheme == ULOGLIN) {
420 G4double a = G4Log(y1);
421 G4double b = (G4Log(y2) - G4Log(y1)) / (x2 - x1);
422 sum += (G4Exp(a) / b) * (G4Exp(b * x2) - G4Exp(b * x1));
423 }
424 else if (aScheme == HISTO || aScheme == CHISTO || aScheme == UHISTO) {
425 sum += y1 * (x2 - x1);
426 }
427 else if (aScheme == LOGLOG || aScheme == CLOGLOG || aScheme == ULOGLOG) {
428 G4double a = G4Log(y1);
429 G4double b = (G4Log(y2) - G4Log(y1)) / (G4Log(x2) - G4Log(x1));
430 sum +=
431 (G4Exp(a) / (b + 1))
432 * (G4Pow::GetInstance()->powA(x2, b + 1) - G4Pow::GetInstance()->powA(x1, b + 1));
433 }
434 else {
435 throw G4HadronicException(
436 __FILE__, __LINE__, "Unknown interpolation scheme in G4ParticleHPVector::Integrate");
437 }
438 }
439 }
440 totalIntegral = sum;
441 }
442
443 inline G4double GetIntegral() // linear interpolation; use with care
444 {
445 if (totalIntegral < -0.5) Integrate();
446 return totalIntegral;
447 }
448
449 inline void SetInterpolationManager(const G4InterpolationManager& aManager)
450 {
451 theManager = aManager;
452 }
453
454 inline const G4InterpolationManager& GetInterpolationManager() const { return theManager; }
455
456 inline void SetInterpolationManager(G4InterpolationManager& aMan) { theManager = aMan; }
457
458 inline void SetScheme(G4int aPoint, const G4InterpolationScheme& aScheme)
459 {
460 theManager.AppendScheme(aPoint, aScheme);
461 }
462
463 inline G4InterpolationScheme GetScheme(G4int anIndex) { return theManager.GetScheme(anIndex); }
464
465 G4double GetMeanX()
466 {
467 G4double result;
468 G4double running = 0;
469 G4double weighted = 0;
470 for (G4int i = 1; i < nEntries; i++) {
471 running +=
472 theInt.GetBinIntegral(theManager.GetScheme(i - 1), theData[i - 1].GetX(),
473 theData[i].GetX(), theData[i - 1].GetY(), theData[i].GetY());
474 weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i - 1),
475 theData[i - 1].GetX(), theData[i].GetX(),
476 theData[i - 1].GetY(), theData[i].GetY());
477 }
478 result = weighted / running;
479 return result;
480 }
481
482 // Finds maximum cross section between two values of kinetic energy
483 G4double GetMaxY(G4double emin, G4double emax);
484
485 /*
486 void Block(G4double aX)
487 {
488 theBlocked.push_back(aX);
489 }
490
491 void Buffer(G4double aX)
492 {
493 theBuffered.push_back(aX);
494 }
495 */
496
497 std::vector<G4double> GetBlocked() { return theBlocked; }
498 std::vector<G4double> GetBuffered() { return theBuffered; }
499
500 // void SetBlocked(const std::vector<G4double> &aBlocked) {theBlocked = aBlocked;}
501 // void SetBuffered(const std::vector<G4double> &aBuffer) {theBuffered = aBuffer;}
502
503 G4double Get15percentBorder();
504 G4double Get50percentBorder();
505
506 private:
507 void Check(G4int i);
508
509 G4bool IsBlocked(G4double aX);
510
511 private:
512 G4ParticleHPInterpolator theLin;
513
514 private:
515 G4double totalIntegral;
516
517 G4ParticleHPDataPoint* theData; // the data
518 G4InterpolationManager theManager; // knows how to interpolate the data.
519 G4double* theIntegral;
520 G4int nEntries;
521 G4int nPoints;
522 G4double label;
523
524 G4ParticleHPInterpolator theInt;
525 G4int Verbose;
526 // debug only
527 G4int isFreed;
528
529 G4ParticleHPHash theHash;
530 G4double maxValue;
531
532 std::vector<G4double> theBlocked;
533 std::vector<G4double> theBuffered;
534 G4double the15percentBorderCash;
535 G4double the50percentBorderCash;
536 };
537
538 #endif
539