| 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 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron transport model.
29 //
30 // 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi)
31 // 080808 bug fix in Sample() and GetXsec() by T. Koi
32 //
33 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
34 //
35 #include "G4ParticleHPVector.hh"
36
37 #include "G4SystemOfUnits.hh"
38 #include "G4Threading.hh"
39
40 // if the ranges do not match, constant extrapolation is used.
41 G4ParticleHPVector& operator+(G4ParticleHPVector& left, G4ParticleHPVector& right)
42 {
43 auto result = new G4ParticleHPVector;
44 G4int j = 0;
45 G4double x;
46 G4double y;
47 G4int running = 0;
48 for (G4int i = 0; i < left.GetVectorLength(); i++) {
49 while (j < right.GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
50 {
51 if (right.GetX(j) < left.GetX(i) * 1.001) {
52 x = right.GetX(j);
53 y = right.GetY(j) + left.GetY(x);
54 result->SetData(running++, x, y);
55 j++;
56 }
57 // else if(std::abs((right.GetX(j)-left.GetX(i))/(left.GetX(i)+right.GetX(j)))>0.001)
58 else if (left.GetX(i) + right.GetX(j) == 0
59 || std::abs((right.GetX(j) - left.GetX(i)) / (left.GetX(i) + right.GetX(j))) > 0.001)
60 {
61 x = left.GetX(i);
62 y = left.GetY(i) + right.GetY(x);
63 result->SetData(running++, x, y);
64 break;
65 }
66 else {
67 break;
68 }
69 }
70 if (j == right.GetVectorLength()) {
71 x = left.GetX(i);
72 y = left.GetY(i) + right.GetY(x);
73 result->SetData(running++, x, y);
74 }
75 }
76 result->ThinOut(0.02);
77 return *result;
78 }
79
80 G4ParticleHPVector::G4ParticleHPVector()
81 {
82 theData = new G4ParticleHPDataPoint[20];
83 nPoints = 20;
84 nEntries = 0;
85 Verbose = 0;
86 theIntegral = nullptr;
87 totalIntegral = -1;
88 isFreed = 0;
89 maxValue = -DBL_MAX;
90 the15percentBorderCash = -DBL_MAX;
91 the50percentBorderCash = -DBL_MAX;
92 label = -DBL_MAX;
93 }
94
95 G4ParticleHPVector::G4ParticleHPVector(G4int n)
96 {
97 nPoints = std::max(n, 20);
98 theData = new G4ParticleHPDataPoint[nPoints];
99 nEntries = 0;
100 Verbose = 0;
101 theIntegral = nullptr;
102 totalIntegral = -1;
103 isFreed = 0;
104 maxValue = -DBL_MAX;
105 the15percentBorderCash = -DBL_MAX;
106 the50percentBorderCash = -DBL_MAX;
107 label = -DBL_MAX;
108 }
109
110 G4ParticleHPVector::~G4ParticleHPVector()
111 {
112 // if(Verbose==1)G4cout <<"G4ParticleHPVector::~G4ParticleHPVector"<<G4endl;
113 delete[] theData;
114 // if(Verbose==1)G4cout <<"Vector: delete theData"<<G4endl;
115 delete[] theIntegral;
116 // if(Verbose==1)G4cout <<"Vector: delete theIntegral"<<G4endl;
117 theHash.Clear();
118 isFreed = 1;
119 }
120
121 G4ParticleHPVector& G4ParticleHPVector::operator=(const G4ParticleHPVector& right)
122 {
123 if (&right == this) return *this;
124
125 G4int i;
126
127 totalIntegral = right.totalIntegral;
128 if (right.theIntegral != nullptr) theIntegral = new G4double[right.nEntries];
129 for (i = 0; i < right.nEntries; i++) {
130 SetPoint(i, right.GetPoint(i)); // copy theData
131 if (right.theIntegral != nullptr) theIntegral[i] = right.theIntegral[i];
132 }
133 theManager = right.theManager;
134 label = right.label;
135
136 Verbose = right.Verbose;
137 the15percentBorderCash = right.the15percentBorderCash;
138 the50percentBorderCash = right.the50percentBorderCash;
139 theHash = right.theHash;
140 return *this;
141 }
142
143 G4double G4ParticleHPVector::GetXsec(G4double e)
144 {
145 if (nEntries == 0) return 0;
146 // if(!theHash.Prepared()) Hash();
147 if (!theHash.Prepared()) {
148 if (G4Threading::IsWorkerThread()) {
149 ;
150 }
151 else {
152 Hash();
153 }
154 }
155 G4int min = theHash.GetMinIndex(e);
156 G4int i;
157 for (i = min; i < nEntries; i++) {
158 // if(theData[i].GetX()>e) break;
159 if (theData[i].GetX() >= e) break;
160 }
161 G4int low = i - 1;
162 G4int high = i;
163 if (i == 0) {
164 low = 0;
165 high = 1;
166 }
167 else if (i == nEntries) {
168 low = nEntries - 2;
169 high = nEntries - 1;
170 }
171 G4double y;
172 if (e < theData[nEntries - 1].GetX()) {
173 // Protect against doubled-up x values
174 // if( (theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001)
175 if (theData[high].GetX() != 0
176 // 080808 TKDB
177 //&&( theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001)
178 && (std::abs((theData[high].GetX() - theData[low].GetX()) / theData[high].GetX())
179 < 0.000001))
180 {
181 y = theData[low].GetY();
182 }
183 else {
184 y = theInt.Interpolate(theManager.GetScheme(high), e, theData[low].GetX(),
185 theData[high].GetX(), theData[low].GetY(), theData[high].GetY());
186 }
187 }
188 else {
189 y = theData[nEntries - 1].GetY();
190 }
191 return y;
192 }
193
194 G4int G4ParticleHPVector::GetEnergyIndex(G4double &e)
195 {
196 // returns energy index of the bin in which ekin is lower then emax and higher than emin for CALENDF
197 // returns energy index of emax for NJOY
198 if (nEntries == 0) return 0;
199 if ( ( ! theHash.Prepared() ) && ( ! G4Threading::IsWorkerThread() ) ) Hash();
200 G4int min = theHash.GetMinIndex(e);
201 G4int i = 0;
202 for (i=min ; i < nEntries; i++) if (theData[i].GetX() >= e) break;
203 return i;
204 }
205
206 void G4ParticleHPVector::Dump()
207 {
208 G4cout << nEntries << G4endl;
209 for (G4int i = 0; i < nEntries; i++) {
210 G4cout << theData[i].GetX() << " ";
211 G4cout << theData[i].GetY() << " ";
212 // if (i!=1&&i==5*(i/5)) G4cout << G4endl;
213 G4cout << G4endl;
214 }
215 G4cout << G4endl;
216 }
217
218 void G4ParticleHPVector::Check(G4int i)
219 {
220 if (i > nEntries)
221 throw G4HadronicException(__FILE__, __LINE__,
222 "Skipped some index numbers in G4ParticleHPVector");
223 if (i == nPoints) {
224 nPoints = static_cast<G4int>(1.2 * nPoints);
225 auto buff = new G4ParticleHPDataPoint[nPoints];
226 for (G4int j = 0; j < nEntries; j++)
227 buff[j] = theData[j];
228 delete[] theData;
229 theData = buff;
230 }
231 if (i == nEntries) nEntries = i + 1;
232 }
233
234 void G4ParticleHPVector::Merge(G4InterpolationScheme aScheme, G4double aValue,
235 G4ParticleHPVector* active, G4ParticleHPVector* passive)
236 {
237 // interpolate between labels according to aScheme, cut at aValue,
238 // continue in unknown areas by substraction of the last difference.
239
240 CleanUp();
241 G4int s_tmp = 0, n = 0, m_tmp = 0;
242 G4ParticleHPVector* tmp;
243 G4int a = s_tmp, p = n, t;
244 while (a < active->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
245 {
246 if (active->GetEnergy(a) <= passive->GetEnergy(p)) {
247 G4double xa = active->GetEnergy(a);
248 G4double yy = theInt.Interpolate(aScheme, aValue, active->GetLabel(), passive->GetLabel(),
249 active->GetXsec(a), passive->GetXsec(xa));
250 SetData(m_tmp, xa, yy);
251 theManager.AppendScheme(m_tmp, active->GetScheme(a));
252 m_tmp++;
253 a++;
254 G4double xp = passive->GetEnergy(p);
255 // if( std::abs(std::abs(xp-xa)/xa)<0.0000001&&a<active->GetVectorLength() )
256 if (xa != 0 && std::abs(std::abs(xp - xa) / xa) < 0.0000001 && a < active->GetVectorLength())
257 {
258 p++;
259 tmp = active;
260 t = a;
261 active = passive;
262 a = p;
263 passive = tmp;
264 p = t;
265 }
266 }
267 else {
268 tmp = active;
269 t = a;
270 active = passive;
271 a = p;
272 passive = tmp;
273 p = t;
274 }
275 }
276
277 G4double deltaX = passive->GetXsec(GetEnergy(m_tmp - 1)) - GetXsec(m_tmp - 1);
278 while (p != passive->GetVectorLength()
279 && passive->GetEnergy(p) <= aValue) // Loop checking, 11.05.2015, T. Koi
280 {
281 G4double anX;
282 anX = passive->GetXsec(p) - deltaX;
283 if (anX > 0) {
284 // if(std::abs(GetEnergy(m-1)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.0000001)
285 if (passive->GetEnergy(p) == 0
286 || std::abs(GetEnergy(m_tmp - 1) - passive->GetEnergy(p)) / passive->GetEnergy(p)
287 > 0.0000001)
288 {
289 SetData(m_tmp, passive->GetEnergy(p), anX);
290 theManager.AppendScheme(m_tmp++, passive->GetScheme(p));
291 }
292 }
293 p++;
294 }
295 // Rebuild the Hash;
296 if (theHash.Prepared()) {
297 ReHash();
298 }
299 }
300
301 void G4ParticleHPVector::ThinOut(G4double precision)
302 {
303 // anything in there?
304 if (GetVectorLength() == 0) return;
305 // make the new vector
306 auto aBuff = new G4ParticleHPDataPoint[nPoints];
307 G4double x, x1, x2, y, y1, y2;
308 G4int count = 0, current = 2, start = 1;
309
310 // First element always goes and is never tested.
311 aBuff[0] = theData[0];
312
313 // Find the rest
314 while (current < GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
315 {
316 x1 = aBuff[count].GetX();
317 y1 = aBuff[count].GetY();
318 x2 = theData[current].GetX();
319 y2 = theData[current].GetY();
320
321 if (x1 - x2 == 0) {
322 // Following block added for avoiding div 0 error on Release + G4FPE_DEBUG
323 for (G4int j = start; j < current; j++) {
324 y = (y2 + y1) / 2.;
325 if (std::abs(y - theData[j].GetY()) > precision * y) {
326 aBuff[++count] = theData[current - 1]; // for this one, everything was fine
327 start = current; // the next candidate
328 break;
329 }
330 }
331 }
332 else {
333 for (G4int j = start; j < current; j++) {
334 x = theData[j].GetX();
335 if (x1 - x2 == 0)
336 y = (y2 + y1) / 2.;
337 else
338 y = theInt.Lin(x, x1, x2, y1, y2);
339 if (std::abs(y - theData[j].GetY()) > precision * y) {
340 aBuff[++count] = theData[current - 1]; // for this one, everything was fine
341 start = current; // the next candidate
342 break;
343 }
344 }
345 }
346 current++;
347 }
348 // The last one also always goes, and is never tested.
349 aBuff[++count] = theData[GetVectorLength() - 1];
350 delete[] theData;
351 theData = aBuff;
352 nEntries = count + 1;
353
354 // Rebuild the Hash;
355 if (theHash.Prepared()) {
356 ReHash();
357 }
358 }
359
360 G4bool G4ParticleHPVector::IsBlocked(G4double aX)
361 {
362 G4bool result = false;
363 std::vector<G4double>::iterator i;
364 for (i = theBlocked.begin(); i != theBlocked.end(); i++) {
365 G4double aBlock = *i;
366 if (std::abs(aX - aBlock) < 0.1 * MeV) {
367 result = true;
368 theBlocked.erase(i);
369 break;
370 }
371 }
372 return result;
373 }
374
375 G4double G4ParticleHPVector::Sample() // Samples X according to distribution Y
376 {
377 G4double result = 0.;
378 G4int j;
379 for (j = 0; j < GetVectorLength(); j++) {
380 if (GetY(j) < 0) SetY(j, 0);
381 }
382
383 if (!theBuffered.empty() && G4UniformRand() < 0.5) {
384 result = theBuffered[0];
385 theBuffered.erase(theBuffered.begin());
386 if (result < GetX(GetVectorLength() - 1)) return result;
387 }
388 if (GetVectorLength() == 1) {
389 result = theData[0].GetX();
390 }
391 else {
392 if (theIntegral == nullptr) {
393 IntegrateAndNormalise();
394 }
395 G4int icounter = 0;
396 G4int icounter_max = 1024;
397 do {
398 icounter++;
399 if (icounter > icounter_max) {
400 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of "
401 << __FILE__ << "." << G4endl;
402 break;
403 }
404 // 080808
405 /*
406 G4double rand;
407 G4double value, test, baseline;
408 baseline = theData[GetVectorLength()-1].GetX()-theData[0].GetX();
409 do
410 {
411 value = baseline*G4UniformRand();
412 value += theData[0].GetX();
413 test = GetY(value)/maxValue;
414 rand = G4UniformRand();
415 }
416 //while(test<rand);
417 while( test < rand && test > 0 );
418 result = value;
419 */
420 G4double rand;
421 G4double value = 0., test;
422 G4int jcounter = 0;
423 G4int jcounter_max = 1024;
424 do {
425 jcounter++;
426 if (jcounter > jcounter_max) {
427 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of "
428 << __FILE__ << "." << G4endl;
429 break;
430 }
431 rand = G4UniformRand();
432 G4int ibin = -1;
433 for (G4int i = 0; i < GetVectorLength(); i++) {
434 if (rand < theIntegral[i]) {
435 ibin = i;
436 break;
437 }
438 }
439 if (ibin < 0) G4cout << "TKDB 080807 " << rand << G4endl;
440 // result
441 rand = G4UniformRand();
442 G4double x1, x2;
443 if (ibin == 0) {
444 x1 = theData[ibin].GetX();
445 value = x1;
446 break;
447 }
448 x1 = theData[ibin - 1].GetX();
449
450 x2 = theData[ibin].GetX();
451 value = rand * (x2 - x1) + x1;
452 //***********************************************************************
453 /*
454 test = GetY ( value ) / std::max ( GetY( ibin-1 ) , GetY ( ibin ) );
455 */
456 //***********************************************************************
457 // EMendoza - Always linear interpolation:
458 G4double y1 = theData[ibin - 1].GetY();
459 G4double y2 = theData[ibin].GetY();
460 G4double mval = (y2 - y1) / (x2 - x1);
461 G4double bval = y1 - mval * x1;
462 test = (mval * value + bval) / std::max(GetY(ibin - 1), GetY(ibin));
463 //***********************************************************************
464 } while (G4UniformRand() > test); // Loop checking, 11.05.2015, T. Koi
465 result = value;
466 // 080807
467 } while (IsBlocked(result)); // Loop checking, 11.05.2015, T. Koi
468 }
469 return result;
470 }
471
472 G4double G4ParticleHPVector::Get15percentBorder()
473 {
474 if (the15percentBorderCash > -DBL_MAX / 2.) return the15percentBorderCash;
475 G4double result;
476 if (GetVectorLength() == 1) {
477 result = theData[0].GetX();
478 the15percentBorderCash = result;
479 }
480 else {
481 if (theIntegral == nullptr) {
482 IntegrateAndNormalise();
483 }
484 G4int i;
485 result = theData[GetVectorLength() - 1].GetX();
486 for (i = 0; i < GetVectorLength(); i++) {
487 if (theIntegral[i] / theIntegral[GetVectorLength() - 1] > 0.15) {
488 result = theData[std::min(i + 1, GetVectorLength() - 1)].GetX();
489 the15percentBorderCash = result;
490 break;
491 }
492 }
493 the15percentBorderCash = result;
494 }
495 return result;
496 }
497
498 G4double G4ParticleHPVector::Get50percentBorder()
499 {
500 if (the50percentBorderCash > -DBL_MAX / 2.) return the50percentBorderCash;
501 G4double result;
502 if (GetVectorLength() == 1) {
503 result = theData[0].GetX();
504 the50percentBorderCash = result;
505 }
506 else {
507 if (theIntegral == nullptr) {
508 IntegrateAndNormalise();
509 }
510 G4int i;
511 G4double x = 0.5;
512 result = theData[GetVectorLength() - 1].GetX();
513 for (i = 0; i < GetVectorLength(); i++) {
514 if (theIntegral[i] / theIntegral[GetVectorLength() - 1] > x) {
515 G4int it;
516 it = i;
517 if (it == GetVectorLength() - 1) {
518 result = theData[GetVectorLength() - 1].GetX();
519 }
520 else {
521 G4double x1, x2, y1, y2;
522 x1 = theIntegral[i - 1] / theIntegral[GetVectorLength() - 1];
523 x2 = theIntegral[i] / theIntegral[GetVectorLength() - 1];
524 y1 = theData[i - 1].GetX();
525 y2 = theData[i].GetX();
526 result = theLin.Lin(x, x1, x2, y1, y2);
527 }
528 the50percentBorderCash = result;
529 break;
530 }
531 }
532 the50percentBorderCash = result;
533 }
534 return result;
535 }
536
537 G4double G4ParticleHPVector::GetMaxY(G4double emin, G4double emax)
538 {
539 G4double xsmax = 0.;
540 if (emin > emax || nEntries == 0) return xsmax;
541 if (emin >= theData[nEntries - 1].GetX()) {
542 xsmax = theData[nEntries - 1].GetY();
543 return xsmax;
544 }
545 if (emax <= theData[0].GetX()) {
546 xsmax = theData[0].GetY();
547 return xsmax;
548 }
549 if (!theHash.Prepared() && !G4Threading::IsWorkerThread()) Hash();
550 // Find the lowest index, low, where x(energy) is higher than emin
551 G4int i = theHash.GetMinIndex(emin);
552 for (; i < nEntries; ++i) {
553 if (theData[i].GetX() >= emin) break;
554 }
555 G4int low = i;
556 // Find the lowest index, high, where x(energy) is higher than emax
557 i = theHash.GetMinIndex(emax);
558 for (; i < nEntries; ++i) {
559 if (theData[i].GetX() >= emax) break;
560 }
561 G4int high = i;
562 xsmax = GetXsec(emin); // Set xsmax as low border
563 // Find the highest cross-section
564 for (i = low; i < high; ++i) {
565 if (xsmax < theData[i].GetY()) xsmax = theData[i].GetY();
566 }
567 // Check if it is smaller than the high border (e.g. as for a monotonic increasing cross-section)
568 G4double highborder = GetXsec(emax);
569 if (xsmax < highborder) xsmax = highborder;
570
571 if (xsmax == 0.) {
572 throw G4HadronicException(__FILE__, __LINE__,
573 "G4ParticleHPVector::GetMaxY : called "
574 "G4Nucleus::GetBiasedThermalNucleus for DBRC, xsmax==0.");
575 }
576
577 return xsmax;
578 }
579