| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans
29 //
30 // 070523 bug fix for G4FPE_DEBUG on by A. How
31 // 080808 bug fix in Sample() and GetXsec() by
32 //
33 // P. Arce, June-2014 Conversion neutron_hp to
34 //
35 #include "G4ParticleHPVector.hh"
36
37 #include "G4SystemOfUnits.hh"
38 #include "G4Threading.hh"
39
40 // if the ranges do not match, constant extrap
41 G4ParticleHPVector& operator+(G4ParticleHPVect
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()
49 while (j < right.GetVectorLength()) // Lo
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.
58 else if (left.GetX(i) + right.GetX(j) ==
59 || std::abs((right.GetX(j) - le
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 <<"G4ParticleHPVe
113 delete[] theData;
114 // if(Verbose==1)G4cout <<"Vector: delete
115 delete[] theIntegral;
116 // if(Verbose==1)G4cout <<"Vector: delete
117 theHash.Clear();
118 isFreed = 1;
119 }
120
121 G4ParticleHPVector& G4ParticleHPVector::operat
122 {
123 if (&right == this) return *this;
124
125 G4int i;
126
127 totalIntegral = right.totalIntegral;
128 if (right.theIntegral != nullptr) theIntegra
129 for (i = 0; i < right.nEntries; i++) {
130 SetPoint(i, right.GetPoint(i)); // copy t
131 if (right.theIntegral != nullptr) theInteg
132 }
133 theManager = right.theManager;
134 label = right.label;
135
136 Verbose = right.Verbose;
137 the15percentBorderCash = right.the15percentB
138 the50percentBorderCash = right.the50percentB
139 theHash = right.theHash;
140 return *this;
141 }
142
143 G4double G4ParticleHPVector::GetXsec(G4double
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].
175 if (theData[high].GetX() != 0
176 // 080808 TKDB
177 //&&( theData[high].GetX()-theData[low
178 && (std::abs((theData[high].GetX() - t
179 < 0.000001))
180 {
181 y = theData[low].GetY();
182 }
183 else {
184 y = theInt.Interpolate(theManager.GetSch
185 theData[high].Get
186 }
187 }
188 else {
189 y = theData[nEntries - 1].GetY();
190 }
191 return y;
192 }
193
194 G4int G4ParticleHPVector::GetEnergyIndex(G4dou
195 {
196 // returns energy index of the bin in which
197 // returns energy index of emax for NJOY
198 if (nEntries == 0) return 0;
199 if ( ( ! theHash.Prepared() ) && ( ! G4Threa
200 G4int min = theHash.GetMinIndex(e);
201 G4int i = 0;
202 for (i=min ; i < nEntries; i++) if (theData[
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 << G4
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 in
223 if (i == nPoints) {
224 nPoints = static_cast<G4int>(1.2 * nPoints
225 auto buff = new G4ParticleHPDataPoint[nPoi
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(G4Interpolation
235 G4ParticleHPVec
236 {
237 // interpolate between labels according to a
238 // continue in unknown areas by substraction
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()) // Lo
245 {
246 if (active->GetEnergy(a) <= passive->GetEn
247 G4double xa = active->GetEnergy(a);
248 G4double yy = theInt.Interpolate(aScheme
249 active-
250 SetData(m_tmp, xa, yy);
251 theManager.AppendScheme(m_tmp, active->G
252 m_tmp++;
253 a++;
254 G4double xp = passive->GetEnergy(p);
255 // if( std::abs(std::abs(xp-xa)/xa)<0.00
256 if (xa != 0 && std::abs(std::abs(xp - xa
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
278 while (p != passive->GetVectorLength()
279 && passive->GetEnergy(p) <= aValue)
280 {
281 G4double anX;
282 anX = passive->GetXsec(p) - deltaX;
283 if (anX > 0) {
284 // if(std::abs(GetEnergy(m-1)-passive->G
285 if (passive->GetEnergy(p) == 0
286 || std::abs(GetEnergy(m_tmp - 1) - p
287 > 0.0000001)
288 {
289 SetData(m_tmp, passive->GetEnergy(p),
290 theManager.AppendScheme(m_tmp++, passi
291 }
292 }
293 p++;
294 }
295 // Rebuild the Hash;
296 if (theHash.Prepared()) {
297 ReHash();
298 }
299 }
300
301 void G4ParticleHPVector::ThinOut(G4double prec
302 {
303 // anything in there?
304 if (GetVectorLength() == 0) return;
305 // make the new vector
306 auto aBuff = new G4ParticleHPDataPoint[nPoin
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 te
311 aBuff[0] = theData[0];
312
313 // Find the rest
314 while (current < GetVectorLength()) // Loop
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 di
323 for (G4int j = start; j < current; j++)
324 y = (y2 + y1) / 2.;
325 if (std::abs(y - theData[j].GetY()) >
326 aBuff[++count] = theData[current - 1
327 start = current; // the next candid
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()) >
340 aBuff[++count] = theData[current - 1
341 start = current; // the next candid
342 break;
343 }
344 }
345 }
346 current++;
347 }
348 // The last one also always goes, and is nev
349 aBuff[++count] = theData[GetVectorLength() -
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
361 {
362 G4bool result = false;
363 std::vector<G4double>::iterator i;
364 for (i = theBlocked.begin(); i != theBlocked
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() // Samp
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()
384 result = theBuffered[0];
385 theBuffered.erase(theBuffered.begin());
386 if (result < GetX(GetVectorLength() - 1))
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 t
401 << __FILE__ << "." << G4endl;
402 break;
403 }
404 // 080808
405 /*
406 G4double rand;
407 G4double value, test, baseline;
408 baseline = theData[GetVectorLeng
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
428 << __FILE__ << "." << G4endl;
429 break;
430 }
431 rand = G4UniformRand();
432 G4int ibin = -1;
433 for (G4int i = 0; i < GetVectorLength(
434 if (rand < theIntegral[i]) {
435 ibin = i;
436 break;
437 }
438 }
439 if (ibin < 0) G4cout << "TKDB 080807 "
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
455 */
456 //************************************
457 // EMendoza - Always linear interpolat
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::ma
463 //************************************
464 } while (G4UniformRand() > test); // Lo
465 result = value;
466 // 080807
467 } while (IsBlocked(result)); // Loop chec
468 }
469 return result;
470 }
471
472 G4double G4ParticleHPVector::Get15percentBorde
473 {
474 if (the15percentBorderCash > -DBL_MAX / 2.)
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].Ge
486 for (i = 0; i < GetVectorLength(); i++) {
487 if (theIntegral[i] / theIntegral[GetVect
488 result = theData[std::min(i + 1, GetVe
489 the15percentBorderCash = result;
490 break;
491 }
492 }
493 the15percentBorderCash = result;
494 }
495 return result;
496 }
497
498 G4double G4ParticleHPVector::Get50percentBorde
499 {
500 if (the50percentBorderCash > -DBL_MAX / 2.)
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].Ge
513 for (i = 0; i < GetVectorLength(); i++) {
514 if (theIntegral[i] / theIntegral[GetVect
515 G4int it;
516 it = i;
517 if (it == GetVectorLength() - 1) {
518 result = theData[GetVectorLength() -
519 }
520 else {
521 G4double x1, x2, y1, y2;
522 x1 = theIntegral[i - 1] / theIntegra
523 x2 = theIntegral[i] / theIntegral[Ge
524 y1 = theData[i - 1].GetX();
525 y2 = theData[i].GetX();
526 result = theLin.Lin(x, x1, x2, y1, y
527 }
528 the50percentBorderCash = result;
529 break;
530 }
531 }
532 the50percentBorderCash = result;
533 }
534 return result;
535 }
536
537 G4double G4ParticleHPVector::GetMaxY(G4double
538 {
539 G4double xsmax = 0.;
540 if (emin > emax || nEntries == 0) return xsm
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::IsW
550 // Find the lowest index, low, where x(energ
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(ener
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
563 // Find the highest cross-section
564 for (i = low; i < high; ++i) {
565 if (xsmax < theData[i].GetY()) xsmax = the
566 }
567 // Check if it is smaller than the high bord
568 G4double highborder = GetXsec(emax);
569 if (xsmax < highborder) xsmax = highborder;
570
571 if (xsmax == 0.) {
572 throw G4HadronicException(__FILE__, __LINE
573 "G4ParticleHPVec
574 "G4Nucleus::GetB
575 }
576
577 return xsmax;
578 }
579