| 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 // 09-May-06 fix in Sample by T. Koi
31 // 080318 Fix Compilation warnings - gcc-4.3.0 by T. Koi
32 // (This fix has a real effect to the code.)
33 // 080409 Fix div0 error with G4FPE by T. Koi
34 // 080612 Fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #1
35 // 080714 Limiting the sum of energy of secondary particles by T. Koi
36 // 080801 Fix div0 error wiht G4FPE and memory leak by T. Koi
37 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
38 //
39 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
40 //
41 // June-2019 - E. Mendoza --> redefinition of the residual mass to consider incident particles
42 // different than neutrons.
43 //
44 // V. Ivanchenko, July-2023 Basic revision of particle HP classes
45 //
46
47 #include "G4ParticleHPContAngularPar.hh"
48
49 #include "G4ParticleDefinition.hh"
50 #include "G4Alpha.hh"
51 #include "G4Deuteron.hh"
52 #include "G4Electron.hh"
53 #include "G4Gamma.hh"
54 #include "G4He3.hh"
55 #include "G4IonTable.hh"
56 #include "G4Neutron.hh"
57 #include "G4NucleiProperties.hh"
58 #include "G4ParticleHPKallbachMannSyst.hh"
59 #include "G4ParticleHPLegendreStore.hh"
60 #include "G4ParticleHPManager.hh"
61 #include "G4ParticleHPVector.hh"
62 #include "G4PhysicalConstants.hh"
63 #include "G4Positron.hh"
64 #include "G4Proton.hh"
65 #include "G4SystemOfUnits.hh"
66 #include "G4Triton.hh"
67
68 #include <set>
69 #include <vector>
70
71 G4ParticleHPContAngularPar::G4ParticleHPContAngularPar(const G4ParticleDefinition* p)
72 {
73 theProjectile = (nullptr == p) ? G4Neutron::Neutron() : p;
74 toBeCached v;
75 fCache.Put(v);
76 if (G4ParticleHPManager::GetInstance()->GetDoNotAdjustFinalState()) adjustResult = false;
77 }
78
79 G4ParticleHPContAngularPar::G4ParticleHPContAngularPar(G4ParticleHPContAngularPar& val)
80 {
81 theEnergy = val.theEnergy;
82 nEnergies = val.nEnergies;
83 nDiscreteEnergies = val.nDiscreteEnergies;
84 nAngularParameters = val.nAngularParameters;
85 theProjectile = val.theProjectile;
86 theManager = val.theManager;
87 theInt = val.theInt;
88 adjustResult = val.adjustResult;
89 theMinEner = val.theMinEner;
90 theMaxEner = val.theMaxEner;
91 theEnergiesTransformed = val.theEnergiesTransformed;
92 theDiscreteEnergies = val.theDiscreteEnergies;
93 theDiscreteEnergiesOwn = val.theDiscreteEnergiesOwn;
94 toBeCached v;
95 fCache.Put(v);
96 const std::size_t esize = nEnergies > 0 ? nEnergies : 1;
97 theAngular = new G4ParticleHPList[esize];
98 for (G4int ie = 0; ie < nEnergies; ++ie) {
99 theAngular[ie].SetLabel(val.theAngular[ie].GetLabel());
100 for (G4int ip = 0; ip < nAngularParameters; ++ip) {
101 theAngular[ie].SetValue(ip, val.theAngular[ie].GetValue(ip));
102 }
103 }
104 }
105
106 G4ParticleHPContAngularPar::~G4ParticleHPContAngularPar()
107 {
108 delete[] theAngular;
109 }
110
111 void G4ParticleHPContAngularPar::Init(std::istream& aDataFile, const G4ParticleDefinition* p)
112 {
113 adjustResult = true;
114 if (G4ParticleHPManager::GetInstance()->GetDoNotAdjustFinalState()) adjustResult = false;
115
116 theProjectile = (nullptr == p) ? G4Neutron::Neutron() : p;
117
118 aDataFile >> theEnergy >> nEnergies >> nDiscreteEnergies >> nAngularParameters;
119 theEnergy *= eV;
120 const std::size_t esize = nEnergies > 0 ? nEnergies : 1;
121 theAngular = new G4ParticleHPList[esize];
122 G4double sEnergy;
123 for (G4int i = 0; i < nEnergies; ++i) {
124 aDataFile >> sEnergy;
125 sEnergy *= eV;
126 theAngular[i].SetLabel(sEnergy);
127 theAngular[i].Init(aDataFile, nAngularParameters, 1.);
128 theMinEner = std::min(theMinEner, sEnergy);
129 theMaxEner = std::max(theMaxEner, sEnergy);
130 }
131 }
132
133 G4ReactionProduct* G4ParticleHPContAngularPar::Sample(G4double anEnergy, G4double massCode,
134 G4double /*targetMass*/, G4int angularRep,
135 G4int /*interpolE*/)
136 {
137 // The following line is needed because it may change between runs by UI command
138 adjustResult = true;
139 if (G4ParticleHPManager::GetInstance()->GetDoNotAdjustFinalState()) adjustResult = false;
140
141 auto result = new G4ReactionProduct;
142 auto Z = static_cast<G4int>(massCode / 1000);
143 auto A = static_cast<G4int>(massCode - 1000 * Z);
144 if (massCode == 0) {
145 result->SetDefinition(G4Gamma::Gamma());
146 }
147 else if (A == 0) {
148 result->SetDefinition(G4Electron::Electron());
149 if (Z == 1) result->SetDefinition(G4Positron::Positron());
150 }
151 else if (A == 1) {
152 result->SetDefinition(G4Neutron::Neutron());
153 if (Z == 1) result->SetDefinition(G4Proton::Proton());
154 }
155 else if (A == 2) {
156 result->SetDefinition(G4Deuteron::Deuteron());
157 }
158 else if (A == 3) {
159 result->SetDefinition(G4Triton::Triton());
160 if (Z == 2) result->SetDefinition(G4He3::He3());
161 }
162 else if (A == 4) {
163 result->SetDefinition(G4Alpha::Alpha());
164 if (Z != 2)
165 throw G4HadronicException(__FILE__, __LINE__,
166 "G4ParticleHPContAngularPar: Unknown ion case 1");
167 }
168 else {
169 result->SetDefinition(G4IonTable::GetIonTable()->GetIon(Z, A, 0));
170 }
171
172 G4int i(0);
173 G4int it(0);
174 G4double fsEnergy(0);
175 G4double cosTh(0);
176 /*
177 G4cout << "G4ParticleHPContAngularPar::Sample E=" << anEnergy <<" Z=" << Z << " A=" << A
178 << " angularRep=" << angularRep << " Nd=" << nDiscreteEnergies
179 << " Ne=" << nEnergies << G4endl;
180 */
181 if (angularRep == 1) {
182 if (nDiscreteEnergies != 0) {
183 // 1st check remaining_energy
184 // if this is the first set it. (How?)
185 if (fCache.Get().fresh) {
186 // Discrete Lines, larger energies come first
187 // Continues Emssions, low to high LAST
188 fCache.Get().remaining_energy =
189 std::max(theAngular[0].GetLabel(), theAngular[nEnergies - 1].GetLabel());
190 fCache.Get().fresh = false;
191 }
192
193 // Cheating for small remaining_energy
194 // Temporary solution
195 if (nDiscreteEnergies == nEnergies) {
196 fCache.Get().remaining_energy =
197 std::max(fCache.Get().remaining_energy,
198 theAngular[nDiscreteEnergies - 1].GetLabel()); // Minimum Line
199 }
200 else {
201 G4double cont_min = 0.0;
202 for (G4int j = nDiscreteEnergies; j < nEnergies; ++j) {
203 cont_min = theAngular[j].GetLabel();
204 if (theAngular[j].GetValue(0) != 0.0) break;
205 }
206 fCache.Get().remaining_energy = std::max(
207 fCache.Get().remaining_energy, std::min(theAngular[nDiscreteEnergies - 1].GetLabel(),
208 cont_min)); // Minimum Line or grid
209 }
210
211 G4double random = G4UniformRand();
212 auto running = new G4double[nEnergies + 1];
213 running[0] = 0.0;
214
215 G4double delta;
216 for (G4int j = 0; j < nDiscreteEnergies; ++j) {
217 delta = 0.0;
218 if (theAngular[j].GetLabel() <= fCache.Get().remaining_energy)
219 delta = theAngular[j].GetValue(0);
220 running[j + 1] = running[j] + delta;
221 }
222
223 G4double tot_prob_DIS = std::max(running[nDiscreteEnergies], 0.0);
224
225 G4double delta1;
226 for (G4int j = nDiscreteEnergies; j < nEnergies; ++j) {
227 delta1 = 0.0;
228 G4double e_low = 0.0;
229 G4double e_high = 0.0;
230 if (theAngular[j].GetLabel() <= fCache.Get().remaining_energy)
231 delta1 = theAngular[j].GetValue(0);
232
233 // To calculate Prob. e_low and e_high should be in eV
234 // There are two cases:
235 // 1: theAngular[nDiscreteEnergies].GetLabel() != 0.0
236 // delta1 should be used between j-1 and j
237 // At j = nDiscreteEnergies (the first) e_low should be set explicitly
238 if (theAngular[j].GetLabel() != 0) {
239 if (j == nDiscreteEnergies) {
240 e_low = 0.0 / eV;
241 }
242 else {
243 if (j < 1) j = 1; // Protection against evaluation of arrays at index j-1
244 e_low = theAngular[j - 1].GetLabel() / eV;
245 }
246 e_high = theAngular[j].GetLabel() / eV;
247 }
248
249 // 2: theAngular[nDiscreteEnergies].GetLabel() == 0.0
250 // delta1 should be used between j and j+1
251 if (theAngular[j].GetLabel() == 0.0) {
252 e_low = theAngular[j].GetLabel() / eV;
253 if (j != nEnergies - 1) {
254 e_high = theAngular[j + 1].GetLabel() / eV;
255 }
256 else {
257 e_high = theAngular[j].GetLabel() / eV;
258 }
259 }
260
261 running[j + 1] = running[j] + ((e_high - e_low) * delta1);
262 }
263 G4double tot_prob_CON = std::max(running[nEnergies] - running[nDiscreteEnergies], 0.0);
264
265 // Give up in the pathological case of null probabilities
266 if (tot_prob_DIS == 0.0 && tot_prob_CON == 0.0) {
267 delete[] running;
268 return result;
269 }
270 // Normalize random
271 random *= (tot_prob_DIS + tot_prob_CON);
272 // 2nd Judge Discrete or not
273
274 // This should be relatively close to 1 For safty
275 if (random <= (tot_prob_DIS / (tot_prob_DIS + tot_prob_CON))
276 || nDiscreteEnergies == nEnergies)
277 {
278 // Discrete Emission
279 for (G4int j = 0; j < nDiscreteEnergies; ++j) {
280 // Here we should use i+1
281 if (random < running[j + 1]) {
282 it = j;
283 break;
284 }
285 }
286 fsEnergy = theAngular[it].GetLabel();
287
288 G4ParticleHPLegendreStore theStore(1);
289 theStore.Init(0, fsEnergy, nAngularParameters);
290 for (G4int j = 0; j < nAngularParameters; ++j) {
291 theStore.SetCoeff(0, j, theAngular[it].GetValue(j));
292 }
293 // use it to sample.
294 cosTh = theStore.SampleMax(fsEnergy);
295 // Done
296 }
297 else {
298 // Continuous emission
299 for (G4int j = nDiscreteEnergies; j < nEnergies; ++j) {
300 // Here we should use i
301 if (random < running[j]) {
302 it = j;
303 break;
304 }
305 }
306
307 if (it < 1) it = 1; // Protection against evaluation of arrays at index it-1
308
309 G4double x1 = running[it - 1];
310 G4double x2 = running[it];
311
312 G4double y1 = 0.0;
313 if (it != nDiscreteEnergies) y1 = theAngular[it - 1].GetLabel();
314 G4double y2 = theAngular[it].GetLabel();
315
316 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it), random, x1, x2, y1, y2);
317
318 G4ParticleHPLegendreStore theStore(2);
319 theStore.Init(0, y1, nAngularParameters);
320 theStore.Init(1, y2, nAngularParameters);
321 theStore.SetManager(theManager);
322 G4int itt;
323 for (G4int j = 0; j < nAngularParameters; ++j) {
324 itt = it;
325 if (it == nDiscreteEnergies) itt = it + 1;
326 // "This case "it-1" has data for Discrete, so we will use an extrpolated values it and
327 // it+1
328 theStore.SetCoeff(0, j, theAngular[itt - 1].GetValue(j));
329 theStore.SetCoeff(1, j, theAngular[itt].GetValue(j));
330 }
331 // use it to sample.
332 cosTh = theStore.SampleMax(fsEnergy);
333
334 // Done
335 }
336
337 // The remaining energy needs to be lowered by the photon energy in *any* case.
338 // Otherwise additional photons with too high energy will be produced - therefore the
339 // adjustResult condition has been removed
340 fCache.Get().remaining_energy -= fsEnergy;
341 delete[] running;
342
343 // end (nDiscreteEnergies != 0) branch
344 }
345 else {
346 // Only continue, TK will clean up
347 if (fCache.Get().fresh) {
348 fCache.Get().remaining_energy = theAngular[nEnergies - 1].GetLabel();
349 fCache.Get().fresh = false;
350 }
351
352 G4double random = G4UniformRand();
353 auto running = new G4double[nEnergies];
354 running[0] = 0;
355 G4double weighted = 0;
356 for (i = 1; i < nEnergies; i++) {
357 running[i] = running[i - 1];
358 if (fCache.Get().remaining_energy >= theAngular[i].GetLabel()) {
359 running[i] += theInt.GetBinIntegral(
360 theManager.GetScheme(i - 1), theAngular[i - 1].GetLabel(), theAngular[i].GetLabel(),
361 theAngular[i - 1].GetValue(0), theAngular[i].GetValue(0));
362 weighted += theInt.GetWeightedBinIntegral(
363 theManager.GetScheme(i - 1), theAngular[i - 1].GetLabel(), theAngular[i].GetLabel(),
364 theAngular[i - 1].GetValue(0), theAngular[i].GetValue(0));
365 }
366 }
367
368 // Cache the mean energy in this distribution
369 if (nEnergies == 1 || running[nEnergies - 1] == 0) {
370 fCache.Get().currentMeanEnergy = 0.0;
371 }
372 else {
373 fCache.Get().currentMeanEnergy = weighted / running[nEnergies - 1];
374 }
375
376 if (nEnergies == 1) it = 0;
377 if (running[nEnergies - 1] != 0) {
378 for (i = 1; i < nEnergies; i++) {
379 it = i;
380 if (random < running[i] / running[nEnergies - 1]) break;
381 }
382 }
383
384 if (running[nEnergies - 1] == 0) it = 0;
385 if (it < nDiscreteEnergies || it == 0) {
386 if (it == 0) {
387 fsEnergy = theAngular[it].GetLabel();
388 G4ParticleHPLegendreStore theStore(1);
389 theStore.Init(0, fsEnergy, nAngularParameters);
390 for (i = 0; i < nAngularParameters; i++) {
391 theStore.SetCoeff(0, i, theAngular[it].GetValue(i));
392 }
393 // use it to sample.
394 cosTh = theStore.SampleMax(fsEnergy);
395 }
396 else {
397 G4double e1, e2;
398 e1 = theAngular[it - 1].GetLabel();
399 e2 = theAngular[it].GetLabel();
400 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it), random,
401 running[it - 1] / running[nEnergies - 1],
402 running[it] / running[nEnergies - 1], e1, e2);
403 // fill a Legendrestore
404 G4ParticleHPLegendreStore theStore(2);
405 theStore.Init(0, e1, nAngularParameters);
406 theStore.Init(1, e2, nAngularParameters);
407 for (i = 0; i < nAngularParameters; i++) {
408 theStore.SetCoeff(0, i, theAngular[it - 1].GetValue(i));
409 theStore.SetCoeff(1, i, theAngular[it].GetValue(i));
410 }
411 // use it to sample.
412 theStore.SetManager(theManager);
413 cosTh = theStore.SampleMax(fsEnergy);
414 }
415 }
416 else { // continuum contribution
417 G4double x1 = running[it - 1] / running[nEnergies - 1];
418 G4double x2 = running[it] / running[nEnergies - 1];
419 G4double y1 = theAngular[it - 1].GetLabel();
420 G4double y2 = theAngular[it].GetLabel();
421 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it), random, x1, x2, y1, y2);
422 G4ParticleHPLegendreStore theStore(2);
423 theStore.Init(0, y1, nAngularParameters);
424 theStore.Init(1, y2, nAngularParameters);
425 theStore.SetManager(theManager);
426 for (i = 0; i < nAngularParameters; i++) {
427 theStore.SetCoeff(0, i, theAngular[it - 1].GetValue(i));
428 theStore.SetCoeff(1, i, theAngular[it].GetValue(i));
429 }
430 // use it to sample.
431 cosTh = theStore.SampleMax(fsEnergy);
432 }
433 delete[] running;
434
435 // The remaining energy needs to be lowered by the photon energy in
436 // *any* case. Otherwise additional photons with too much energy will be
437 // produced - therefore the adjustResult condition has been removed
438
439 fCache.Get().remaining_energy -= fsEnergy;
440 // end if (nDiscreteEnergies != 0)
441 }
442 // end of (angularRep == 1) branch
443 }
444 else if (angularRep == 2) {
445 // first get the energy (already the right for this incoming energy)
446 G4int j;
447 auto running = new G4double[nEnergies];
448 running[0] = 0;
449 G4double weighted = 0;
450 for (j = 1; j < nEnergies; ++j) {
451 if (j != 0) running[j] = running[j - 1];
452 running[j] += theInt.GetBinIntegral(theManager.GetScheme(j - 1), theAngular[j - 1].GetLabel(),
453 theAngular[j].GetLabel(), theAngular[j - 1].GetValue(0),
454 theAngular[j].GetValue(0));
455 weighted += theInt.GetWeightedBinIntegral(
456 theManager.GetScheme(j - 1), theAngular[j - 1].GetLabel(), theAngular[j].GetLabel(),
457 theAngular[j - 1].GetValue(0), theAngular[j].GetValue(0));
458 }
459
460 // Cache the mean energy in this distribution
461 if (nEnergies == 1)
462 fCache.Get().currentMeanEnergy = 0.0;
463 else
464 fCache.Get().currentMeanEnergy = weighted / running[nEnergies - 1];
465
466 G4int itt(0);
467 G4double randkal = G4UniformRand();
468 for (j = 1; j < nEnergies; ++j) {
469 itt = j;
470 if (randkal*running[nEnergies - 1] < running[j]) break;
471 }
472
473 // Interpolate the secondary energy
474 G4double x, x1, x2, y1, y2;
475 if (itt == 0) itt = 1;
476 x = randkal * running[nEnergies - 1];
477 x1 = running[itt - 1];
478 x2 = running[itt];
479 G4double compoundFraction;
480 // interpolate energy
481 y1 = theAngular[itt - 1].GetLabel();
482 y2 = theAngular[itt].GetLabel();
483 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(itt - 1), x, x1, x2, y1, y2);
484
485 // For theta, interpolate the compoundFractions
486 G4double cLow = theAngular[itt - 1].GetValue(1);
487 G4double cHigh = theAngular[itt].GetValue(1);
488 compoundFraction = theInt.Interpolate(theManager.GetScheme(itt), fsEnergy, y1, y2, cLow, cHigh);
489
490 if (compoundFraction > 1.0)
491 compoundFraction = 1.0; // Protection against unphysical interpolation
492
493 delete[] running;
494
495 // get cosTh
496 G4double incidentEnergy = anEnergy;
497 G4double incidentMass = theProjectile->GetPDGMass();
498 G4double productEnergy = fsEnergy;
499 G4double productMass = result->GetMass();
500 auto targetZ = G4int(fCache.Get().theTargetCode / 1000);
501 auto targetA = G4int(fCache.Get().theTargetCode - 1000 * targetZ);
502
503 // To correspond to natural composition (-nat-) data files.
504 if (targetA == 0) targetA = G4int(fCache.Get().theTarget->GetMass() / amu_c2 + 0.5);
505 G4double targetMass = fCache.Get().theTarget->GetMass();
506 auto incidentA = G4int(incidentMass / amu_c2 + 0.5);
507 auto incidentZ = G4int(theProjectile->GetPDGCharge() + 0.5);
508 G4int residualA = targetA + incidentA - A;
509 G4int residualZ = targetZ + incidentZ - Z;
510 G4double residualMass = G4NucleiProperties::GetNuclearMass(residualA, residualZ);
511
512 G4ParticleHPKallbachMannSyst theKallbach(
513 compoundFraction, incidentEnergy, incidentMass, productEnergy, productMass, residualMass,
514 residualA, residualZ, targetMass, targetA, targetZ, incidentA, incidentZ, A, Z);
515 cosTh = theKallbach.Sample(anEnergy);
516 // end (angularRep == 2) branch
517 }
518 else if (angularRep > 10 && angularRep < 16) {
519 G4double random = G4UniformRand();
520 auto running = new G4double[nEnergies];
521 running[0] = 0;
522 G4double weighted = 0;
523 for (i = 1; i < nEnergies; ++i) {
524 if (i != 0) running[i] = running[i - 1];
525 running[i] += theInt.GetBinIntegral(theManager.GetScheme(i - 1), theAngular[i - 1].GetLabel(),
526 theAngular[i].GetLabel(), theAngular[i - 1].GetValue(0),
527 theAngular[i].GetValue(0));
528 weighted += theInt.GetWeightedBinIntegral(
529 theManager.GetScheme(i - 1), theAngular[i - 1].GetLabel(), theAngular[i].GetLabel(),
530 theAngular[i - 1].GetValue(0), theAngular[i].GetValue(0));
531 }
532
533 // Cache the mean energy in this distribution
534 if (nEnergies == 1)
535 fCache.Get().currentMeanEnergy = 0.0;
536 else
537 fCache.Get().currentMeanEnergy = weighted / running[nEnergies - 1];
538
539 if (nEnergies == 1) it = 0;
540 for (i = 1; i < nEnergies; i++) {
541 it = i;
542 if (random < running[i] / running[nEnergies - 1]) break;
543 }
544
545 if (it < nDiscreteEnergies || it == 0) {
546 if (it == 0) {
547 fsEnergy = theAngular[0].GetLabel();
548 G4ParticleHPVector theStore;
549 G4int aCounter = 0;
550 for (G4int j = 1; j < nAngularParameters; j += 2) {
551 theStore.SetX(aCounter, theAngular[0].GetValue(j));
552 theStore.SetY(aCounter, theAngular[0].GetValue(j + 1));
553 aCounter++;
554 }
555 G4InterpolationManager aMan;
556 aMan.Init(angularRep - 10, nAngularParameters - 1);
557 theStore.SetInterpolationManager(aMan);
558 cosTh = theStore.Sample();
559 }
560 else {
561 fsEnergy = theAngular[it].GetLabel();
562 G4ParticleHPVector theStore;
563 G4InterpolationManager aMan;
564 aMan.Init(angularRep - 10, nAngularParameters - 1);
565 theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
566 G4InterpolationScheme currentScheme = theManager.GetInverseScheme(it);
567 G4int aCounter = 0;
568 for (G4int j = 1; j < nAngularParameters; j += 2) {
569 theStore.SetX(aCounter, theAngular[it].GetValue(j));
570 theStore.SetY(aCounter, theInt.Interpolate(currentScheme, random,
571 running[it - 1] / running[nEnergies - 1],
572 running[it] / running[nEnergies - 1],
573 theAngular[it - 1].GetValue(j + 1),
574 theAngular[it].GetValue(j + 1)));
575 ++aCounter;
576 }
577 cosTh = theStore.Sample();
578 }
579 }
580 else {
581 G4double x1 = running[it - 1] / running[nEnergies - 1];
582 G4double x2 = running[it] / running[nEnergies - 1];
583 G4double y1 = theAngular[it - 1].GetLabel();
584 G4double y2 = theAngular[it].GetLabel();
585 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it), random, x1, x2, y1, y2);
586 G4ParticleHPVector theBuff1;
587 G4ParticleHPVector theBuff2;
588 G4InterpolationManager aMan;
589 aMan.Init(angularRep - 10, nAngularParameters - 1);
590
591 G4int j;
592 for (i = 0, j = 1; i < nAngularParameters; i++, j += 2) {
593 theBuff1.SetX(i, theAngular[it - 1].GetValue(j));
594 theBuff1.SetY(i, theAngular[it - 1].GetValue(j + 1));
595 theBuff2.SetX(i, theAngular[it].GetValue(j));
596 theBuff2.SetY(i, theAngular[it].GetValue(j + 1));
597 }
598
599 G4ParticleHPVector theStore;
600 theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
601 x1 = y1;
602 x2 = y2;
603 G4double x, y;
604 for (i = 0; i < theBuff1.GetVectorLength(); i++) {
605 x = theBuff1.GetX(i); // costh binning identical
606 y1 = theBuff1.GetY(i);
607 y2 = theBuff2.GetY(i);
608 y = theInt.Interpolate(theManager.GetScheme(it), fsEnergy, theAngular[it - 1].GetLabel(),
609 theAngular[it].GetLabel(), y1, y2);
610 theStore.SetX(i, x);
611 theStore.SetY(i, y);
612 }
613 cosTh = theStore.Sample();
614 }
615 delete[] running;
616 }
617 else {
618 throw G4HadronicException(__FILE__, __LINE__,
619 "G4ParticleHPContAngularPar::Sample: Unknown angular representation");
620 }
621 //G4cout << " Efin=" << fsEnergy << G4endl;
622 result->SetKineticEnergy(fsEnergy);
623
624 G4double phi = twopi * G4UniformRand();
625 if(cosTh > 1.0) { cosTh = 1.0; }
626 else if (cosTh < -1.0) { cosTh = -1.0; }
627 G4double sinth = std::sqrt((1.0 - cosTh)*(1.0 + cosTh));
628 G4double mtot = result->GetTotalMomentum();
629 G4ThreeVector tempVector(mtot * sinth * std::cos(phi), mtot * sinth * std::sin(phi), mtot * cosTh);
630 result->SetMomentum(tempVector);
631 return result;
632 }
633
634 void G4ParticleHPContAngularPar::PrepareTableInterpolation()
635 {
636 // Discrete energies: store own energies in a map for faster searching
637 //
638 // The data files sometimes have identical discrete energies (likely typos)
639 // which would lead to overwriting the already existing index and hence
640 // creating a hole in the lookup table.
641 // No attempt is made here to correct for the energies - rather an epsilon
642 // is subtracted from the energy in order to uniquely identify the line
643
644 for (G4int ie = 0; ie < nDiscreteEnergies; ie++) {
645 // check if energy is already present and subtract epsilon if that's the case
646 G4double myE = theAngular[ie].GetLabel();
647 while (theDiscreteEnergiesOwn.find(myE) != theDiscreteEnergiesOwn.end()) {
648 myE -= 1e-6;
649 }
650 theDiscreteEnergiesOwn[myE] = ie;
651 }
652 return;
653 }
654
655 void G4ParticleHPContAngularPar::BuildByInterpolation(G4double anEnergy,
656 G4InterpolationScheme aScheme,
657 G4ParticleHPContAngularPar& angpar1,
658 G4ParticleHPContAngularPar& angpar2)
659 {
660 G4int ie, ie1, ie2, ie1Prev, ie2Prev;
661 // Only rebuild the interpolation table if there is a new interaction.
662 // For several subsequent samplings of final state particles in the same
663 // interaction the existing table should be used
664 if (!fCache.Get().fresh) return;
665
666 // Make copies of angpar1 and angpar2. Since these are given by reference
667 // it can not be excluded that one of them is "this". Hence this code uses
668 // potentially the old "this" for creating the new this - which leads to
669 // memory corruption if the old is not stored as separarte object for lookup
670 const G4ParticleHPContAngularPar copyAngpar1(angpar1), copyAngpar2(angpar2);
671
672 nAngularParameters = copyAngpar1.nAngularParameters;
673 theManager = copyAngpar1.theManager;
674 theEnergy = anEnergy;
675 theMinEner = DBL_MAX; // min and max will be re-calculated after interpolation
676 theMaxEner = -DBL_MAX;
677
678 // The two discrete sets must be merged. A vector holds the temporary data to
679 // be copied to the array in the end. Since the G4ParticleHPList class
680 // contains pointers, can't simply assign elements of this type. Each member
681 // needs to call the explicit Set() method instead.
682
683 // First, average probabilities for those lines that are in both sets
684 const std::map<G4double, G4int> discEnerOwn1 = copyAngpar1.GetDiscreteEnergiesOwn();
685 const std::map<G4double, G4int> discEnerOwn2 = copyAngpar2.GetDiscreteEnergiesOwn();
686 std::map<G4double, G4int>::const_iterator itedeo1;
687 std::map<G4double, G4int>::const_iterator itedeo2;
688 std::vector<G4ParticleHPList*> vAngular(discEnerOwn1.size());
689 G4double discEner1;
690 for (itedeo1 = discEnerOwn1.cbegin(); itedeo1 != discEnerOwn1.cend(); ++itedeo1) {
691 discEner1 = itedeo1->first;
692 if (discEner1 < theMinEner) {
693 theMinEner = discEner1;
694 }
695 if (discEner1 > theMaxEner) {
696 theMaxEner = discEner1;
697 }
698 ie1 = itedeo1->second;
699 itedeo2 = discEnerOwn2.find(discEner1);
700 if (itedeo2 == discEnerOwn2.cend()) {
701 ie2 = -1;
702 }
703 else {
704 ie2 = itedeo2->second;
705 }
706 vAngular[ie1] = new G4ParticleHPList();
707 vAngular[ie1]->SetLabel(copyAngpar1.theAngular[ie1].GetLabel());
708 G4double val1, val2;
709 for (G4int ip = 0; ip < nAngularParameters; ++ip) {
710 val1 = copyAngpar1.theAngular[ie1].GetValue(ip);
711 if (ie2 != -1) {
712 val2 = copyAngpar2.theAngular[ie2].GetValue(ip);
713 }
714 else {
715 val2 = 0.;
716 }
717 G4double value = theInt.Interpolate(aScheme, anEnergy, copyAngpar1.theEnergy,
718 copyAngpar2.theEnergy, val1, val2);
719 vAngular[ie1]->SetValue(ip, value);
720 }
721 } // itedeo1 loop
722
723 // Add the ones in set2 but not in set1
724 std::vector<G4ParticleHPList*>::const_iterator itv;
725 G4double discEner2;
726 for (itedeo2 = discEnerOwn2.cbegin(); itedeo2 != discEnerOwn2.cend(); ++itedeo2) {
727 discEner2 = itedeo2->first;
728 ie2 = itedeo2->second;
729 G4bool notFound = true;
730 itedeo1 = discEnerOwn1.find(discEner2);
731 if (itedeo1 != discEnerOwn1.cend()) {
732 notFound = false;
733 }
734 if (notFound) {
735 // not yet in list
736 if (discEner2 < theMinEner) {
737 theMinEner = discEner2;
738 }
739 if (discEner2 > theMaxEner) {
740 theMaxEner = discEner2;
741 }
742 // find position to insert
743 G4bool isInserted = false;
744 ie = 0;
745 for (itv = vAngular.cbegin(); itv != vAngular.cend(); ++itv, ++ie) {
746 if (discEner2 > (*itv)->GetLabel()) {
747 itv = vAngular.insert(itv, new G4ParticleHPList);
748 (*itv)->SetLabel(copyAngpar2.theAngular[ie2].GetLabel());
749 isInserted = true;
750 break;
751 }
752 }
753 if (!isInserted) {
754 ie = (G4int)vAngular.size();
755 vAngular.push_back(new G4ParticleHPList);
756 vAngular[ie]->SetLabel(copyAngpar2.theAngular[ie2].GetLabel());
757 isInserted = true;
758 }
759
760 G4double val1, val2;
761 for (G4int ip = 0; ip < nAngularParameters; ++ip) {
762 val1 = 0;
763 val2 = copyAngpar2.theAngular[ie2].GetValue(ip);
764 G4double value = theInt.Interpolate(aScheme, anEnergy, copyAngpar1.theEnergy,
765 copyAngpar2.theEnergy, val1, val2);
766 vAngular[ie]->SetValue(ip, value);
767 }
768 } // end if(notFound)
769 } // end loop on itedeo2
770
771 // Store new discrete list
772 nDiscreteEnergies = (G4int)vAngular.size();
773 delete[] theAngular;
774 theAngular = nullptr;
775 if (nDiscreteEnergies > 0) {
776 theAngular = new G4ParticleHPList[nDiscreteEnergies];
777 }
778 theDiscreteEnergiesOwn.clear();
779 theDiscreteEnergies.clear();
780 for (ie = 0; ie < nDiscreteEnergies; ++ie) {
781 theAngular[ie].SetLabel(vAngular[ie]->GetLabel());
782 for (G4int ip = 0; ip < nAngularParameters; ++ip) {
783 theAngular[ie].SetValue(ip, vAngular[ie]->GetValue(ip));
784 }
785 theDiscreteEnergiesOwn[theAngular[ie].GetLabel()] = ie;
786 theDiscreteEnergies.insert(theAngular[ie].GetLabel());
787 }
788
789 // The continuous energies need to be made from scratch like the discrete
790 // ones. Therefore the re-assignemnt of theAngular needs to be done
791 // after the continuous energy set is also finalized. Only then the
792 // total number of nEnergies is known and the array can be allocated.
793
794 // Get minimum and maximum energy interpolating
795 // Don't use theMinEner or theMaxEner here, since the transformed energies
796 // need the interpolated range from the original Angpar
797 G4double interMinEner = copyAngpar1.GetMinEner()
798 + (theEnergy - copyAngpar1.GetEnergy())
799 * (copyAngpar2.GetMinEner() - copyAngpar1.GetMinEner())
800 / (copyAngpar2.GetEnergy() - copyAngpar1.GetEnergy());
801 G4double interMaxEner = copyAngpar1.GetMaxEner()
802 + (theEnergy - copyAngpar1.GetEnergy())
803 * (copyAngpar2.GetMaxEner() - copyAngpar1.GetMaxEner())
804 / (copyAngpar2.GetEnergy() - copyAngpar1.GetEnergy());
805
806 // Loop to energies of new set
807 theEnergiesTransformed.clear();
808
809 G4int nEnergies1 = copyAngpar1.GetNEnergies();
810 G4int nDiscreteEnergies1 = copyAngpar1.GetNDiscreteEnergies();
811 G4double minEner1 = copyAngpar1.GetMinEner();
812 G4double maxEner1 = copyAngpar1.GetMaxEner();
813 G4int nEnergies2 = copyAngpar2.GetNEnergies();
814 G4int nDiscreteEnergies2 = copyAngpar2.GetNDiscreteEnergies();
815 G4double minEner2 = copyAngpar2.GetMinEner();
816 G4double maxEner2 = copyAngpar2.GetMaxEner();
817
818 // First build the list of transformed energies normalized
819 // to the new min max by assuming that the min-max range of
820 // each set would be scalable to the new, interpolated min
821 // max range
822
823 G4double e1(0.);
824 G4double eTNorm1(0.);
825 for (ie1 = nDiscreteEnergies1; ie1 < nEnergies1; ++ie1) {
826 e1 = copyAngpar1.theAngular[ie1].GetLabel();
827 eTNorm1 = (e1 - minEner1);
828 if (maxEner1 != minEner1) eTNorm1 /= (maxEner1 - minEner1);
829 if (eTNorm1 >= 0 && eTNorm1 <= 1) theEnergiesTransformed.insert(eTNorm1);
830 }
831
832 G4double e2(0.);
833 G4double eTNorm2(0.);
834 for (ie2 = nDiscreteEnergies2; ie2 < nEnergies2; ++ie2) {
835 e2 = copyAngpar2.theAngular[ie2].GetLabel();
836 eTNorm2 = (e2 - minEner2);
837 if (maxEner2 != minEner2) eTNorm2 /= (maxEner2 - minEner2);
838 if (eTNorm2 >= 0 && eTNorm2 <= 1) theEnergiesTransformed.insert(eTNorm2);
839 }
840
841 // Now the list of energies is complete
842 nEnergies = nDiscreteEnergies + (G4int)theEnergiesTransformed.size();
843
844 // Create final array of angular parameters
845 const std::size_t esize = nEnergies > 0 ? nEnergies : 1;
846 auto theNewAngular = new G4ParticleHPList[esize];
847
848 // Copy discrete energies and interpolated parameters to new array
849
850 if (theAngular != nullptr) {
851 for (ie = 0; ie < nDiscreteEnergies; ++ie) {
852 theNewAngular[ie].SetLabel(theAngular[ie].GetLabel());
853 for (G4int ip = 0; ip < nAngularParameters; ++ip) {
854 theNewAngular[ie].SetValue(ip, theAngular[ie].GetValue(ip));
855 }
856 }
857 delete[] theAngular;
858 }
859 theAngular = theNewAngular;
860
861 // Interpolate the continuous energies for new array
862 auto iteet = theEnergiesTransformed.begin();
863
864 G4double e1Interp(0.);
865 G4double e2Interp(0.);
866 for (ie = nDiscreteEnergies; ie < nEnergies; ++ie, ++iteet) {
867 G4double eT = (*iteet);
868
869 //--- Use eT1 = eT: Get energy and parameters of copyAngpar1 for this eT
870 e1Interp = (maxEner1 - minEner1) * eT + minEner1;
871 //----- Get parameter value corresponding to this e1Interp
872 for (ie1 = nDiscreteEnergies1; ie1 < nEnergies1; ++ie1) {
873 if ((copyAngpar1.theAngular[ie1].GetLabel() - e1Interp) > 1.E-10 * e1Interp) break;
874 }
875 ie1Prev = ie1 - 1;
876 if (ie1 == 0) ++ie1Prev;
877 if (ie1 == nEnergies1) {
878 ie1--;
879 ie1Prev = ie1;
880 }
881
882 //--- Use eT2 = eT: Get energy and parameters of copyAngpar2 for this eT
883 e2Interp = (maxEner2 - minEner2) * eT + minEner2;
884 //----- Get parameter value corresponding to this e2Interp
885 for (ie2 = nDiscreteEnergies2; ie2 < nEnergies2; ++ie2) {
886 if ((copyAngpar2.theAngular[ie2].GetLabel() - e2Interp) > 1.E-10 * e2Interp) break;
887 }
888 ie2Prev = ie2 - 1;
889 if (ie2 == 0) ++ie2Prev;
890 if (ie2 == nEnergies2) {
891 ie2--;
892 ie2Prev = ie2;
893 }
894
895 //---- Energy corresponding to energy transformed
896 G4double eN = (interMaxEner - interMinEner) * eT + interMinEner;
897
898 theAngular[ie].SetLabel(eN);
899 if (eN < theMinEner) {
900 theMinEner = eN;
901 }
902 if (eN > theMaxEner) {
903 theMaxEner = eN;
904 }
905
906 G4double val1(0.);
907 G4double val2(0.);
908 G4double value(0.);
909 for (G4int ip = 0; ip < nAngularParameters; ++ip) {
910 val1 = theInt.Interpolate2(
911 theManager.GetScheme(ie), e1Interp, copyAngpar1.theAngular[ie1Prev].GetLabel(),
912 copyAngpar1.theAngular[ie1].GetLabel(), copyAngpar1.theAngular[ie1Prev].GetValue(ip),
913 copyAngpar1.theAngular[ie1].GetValue(ip))
914 * (maxEner1 - minEner1);
915 val2 = theInt.Interpolate2(
916 theManager.GetScheme(ie), e2Interp, copyAngpar2.theAngular[ie2Prev].GetLabel(),
917 copyAngpar2.theAngular[ie2].GetLabel(), copyAngpar2.theAngular[ie2Prev].GetValue(ip),
918 copyAngpar2.theAngular[ie2].GetValue(ip))
919 * (maxEner2 - minEner2);
920
921 value = theInt.Interpolate(aScheme, anEnergy, copyAngpar1.theEnergy, copyAngpar2.theEnergy,
922 val1, val2);
923 if (interMaxEner != interMinEner) {
924 value /= (interMaxEner - interMinEner);
925 }
926 else if (value != 0) {
927 throw G4HadronicException(__FILE__, __LINE__,
928 "G4ParticleHPContAngularPar::PrepareTableInterpolation "
929 "interMaxEner == interMinEner and value != 0.");
930 }
931 theAngular[ie].SetValue(ip, value);
932 }
933 } // end loop on nDiscreteEnergies
934
935 for (itv = vAngular.cbegin(); itv != vAngular.cend(); ++itv)
936 delete (*itv);
937 }
938
939 void G4ParticleHPContAngularPar::Dump() const
940 {
941 G4cout << theEnergy << " " << nEnergies << " " << nDiscreteEnergies << " " << nAngularParameters
942 << G4endl;
943
944 for (G4int ii = 0; ii < nEnergies; ++ii)
945 theAngular[ii].Dump();
946 }
947