| 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 // 09-May-06 fix in Sample by T. Koi
31 // 080318 Fix Compilation warnings - gcc-4.3.0
32 // (This fix has a real effect to the c
33 // 080409 Fix div0 error with G4FPE by T. Koi
34 // 080612 Fix contribution from Benoit Pirard
35 // 080714 Limiting the sum of energy of second
36 // 080801 Fix div0 error wiht G4FPE and memory
37 // 081024 G4NucleiPropertiesTable:: to G4Nucle
38 //
39 // P. Arce, June-2014 Conversion neutron_hp to
40 //
41 // June-2019 - E. Mendoza --> redefinition of
42 // different than neutrons.
43 //
44 // V. Ivanchenko, July-2023 Basic revision of
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::G4ParticleHPContAn
72 {
73 theProjectile = (nullptr == p) ? G4Neutron::
74 toBeCached v;
75 fCache.Put(v);
76 if (G4ParticleHPManager::GetInstance()->GetD
77 }
78
79 G4ParticleHPContAngularPar::G4ParticleHPContAn
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.theEnergiesTran
92 theDiscreteEnergies = val.theDiscreteEnergie
93 theDiscreteEnergiesOwn = val.theDiscreteEner
94 toBeCached v;
95 fCache.Put(v);
96 const std::size_t esize = nEnergies > 0 ? nE
97 theAngular = new G4ParticleHPList[esize];
98 for (G4int ie = 0; ie < nEnergies; ++ie) {
99 theAngular[ie].SetLabel(val.theAngular[ie]
100 for (G4int ip = 0; ip < nAngularParameters
101 theAngular[ie].SetValue(ip, val.theAngul
102 }
103 }
104 }
105
106 G4ParticleHPContAngularPar::~G4ParticleHPContA
107 {
108 delete[] theAngular;
109 }
110
111 void G4ParticleHPContAngularPar::Init(std::ist
112 {
113 adjustResult = true;
114 if (G4ParticleHPManager::GetInstance()->GetD
115
116 theProjectile = (nullptr == p) ? G4Neutron::
117
118 aDataFile >> theEnergy >> nEnergies >> nDisc
119 theEnergy *= eV;
120 const std::size_t esize = nEnergies > 0 ? nE
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, nAngularPara
128 theMinEner = std::min(theMinEner, sEnergy)
129 theMaxEner = std::max(theMaxEner, sEnergy)
130 }
131 }
132
133 G4ReactionProduct* G4ParticleHPContAngularPar:
134
135
136 {
137 // The following line is needed because it m
138 adjustResult = true;
139 if (G4ParticleHPManager::GetInstance()->GetD
140
141 auto result = new G4ReactionProduct;
142 auto Z = static_cast<G4int>(massCode / 1000)
143 auto A = static_cast<G4int>(massCode - 1000
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(G4Positr
150 }
151 else if (A == 1) {
152 result->SetDefinition(G4Neutron::Neutron()
153 if (Z == 1) result->SetDefinition(G4Proton
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::H
161 }
162 else if (A == 4) {
163 result->SetDefinition(G4Alpha::Alpha());
164 if (Z != 2)
165 throw G4HadronicException(__FILE__, __LI
166 "G4ParticleHPC
167 }
168 else {
169 result->SetDefinition(G4IonTable::GetIonTa
170 }
171
172 G4int i(0);
173 G4int it(0);
174 G4double fsEnergy(0);
175 G4double cosTh(0);
176 /*
177 G4cout << "G4ParticleHPContAngularPar::Sampl
178 << " angularRep=" << angularRep << "
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 com
187 // Continues Emssions, low to high
188 fCache.Get().remaining_energy =
189 std::max(theAngular[0].GetLabel(), t
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_ener
198 theAngular[nDiscreteEnergie
199 }
200 else {
201 G4double cont_min = 0.0;
202 for (G4int j = nDiscreteEnergies; j <
203 cont_min = theAngular[j].GetLabel();
204 if (theAngular[j].GetValue(0) != 0.0
205 }
206 fCache.Get().remaining_energy = std::m
207 fCache.Get().remaining_energy, std::
208
209 }
210
211 G4double random = G4UniformRand();
212 auto running = new G4double[nEnergies +
213 running[0] = 0.0;
214
215 G4double delta;
216 for (G4int j = 0; j < nDiscreteEnergies;
217 delta = 0.0;
218 if (theAngular[j].GetLabel() <= fCache
219 delta = theAngular[j].GetValue(0);
220 running[j + 1] = running[j] + delta;
221 }
222
223 G4double tot_prob_DIS = std::max(running
224
225 G4double delta1;
226 for (G4int j = nDiscreteEnergies; j < nE
227 delta1 = 0.0;
228 G4double e_low = 0.0;
229 G4double e_high = 0.0;
230 if (theAngular[j].GetLabel() <= fCache
231 delta1 = theAngular[j].GetValue(0);
232
233 // To calculate Prob. e_low and e_high
234 // There are two cases:
235 // 1: theAngular[nDiscreteEnergies].Ge
236 // delta1 should be used between j-
237 // At j = nDiscreteEnergies (the fi
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 a
244 e_low = theAngular[j - 1].GetLabel
245 }
246 e_high = theAngular[j].GetLabel() /
247 }
248
249 // 2: theAngular[nDiscreteEnergies].Ge
250 // delta1 should be used between j
251 if (theAngular[j].GetLabel() == 0.0) {
252 e_low = theAngular[j].GetLabel() / e
253 if (j != nEnergies - 1) {
254 e_high = theAngular[j + 1].GetLabe
255 }
256 else {
257 e_high = theAngular[j].GetLabel()
258 }
259 }
260
261 running[j + 1] = running[j] + ((e_high
262 }
263 G4double tot_prob_CON = std::max(running
264
265 // Give up in the pathological case of n
266 if (tot_prob_DIS == 0.0 && tot_prob_CON
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
275 if (random <= (tot_prob_DIS / (tot_prob_
276 || nDiscreteEnergies == nEnergies)
277 {
278 // Discrete Emission
279 for (G4int j = 0; j < nDiscreteEnergie
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, nAngularPar
290 for (G4int j = 0; j < nAngularParamete
291 theStore.SetCoeff(0, j, theAngular[i
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 <
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 aga
308
309 G4double x1 = running[it - 1];
310 G4double x2 = running[it];
311
312 G4double y1 = 0.0;
313 if (it != nDiscreteEnergies) y1 = theA
314 G4double y2 = theAngular[it].GetLabel(
315
316 fsEnergy = theInt.Interpolate(theManag
317
318 G4ParticleHPLegendreStore theStore(2);
319 theStore.Init(0, y1, nAngularParameter
320 theStore.Init(1, y2, nAngularParameter
321 theStore.SetManager(theManager);
322 G4int itt;
323 for (G4int j = 0; j < nAngularParamete
324 itt = it;
325 if (it == nDiscreteEnergies) itt = i
326 // "This case "it-1" has data for Di
327 // it+1
328 theStore.SetCoeff(0, j, theAngular[i
329 theStore.SetCoeff(1, j, theAngular[i
330 }
331 // use it to sample.
332 cosTh = theStore.SampleMax(fsEnergy);
333
334 // Done
335 }
336
337 // The remaining energy needs to be lowe
338 // Otherwise additional photons with too
339 // adjustResult condition has been remov
340 fCache.Get().remaining_energy -= fsEnerg
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 = theAng
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 >= t
359 running[i] += theInt.GetBinIntegral(
360 theManager.GetScheme(i - 1), theAn
361 theAngular[i - 1].GetValue(0), the
362 weighted += theInt.GetWeightedBinInt
363 theManager.GetScheme(i - 1), theAn
364 theAngular[i - 1].GetValue(0), the
365 }
366 }
367
368 // Cache the mean energy in this distrib
369 if (nEnergies == 1 || running[nEnergies
370 fCache.Get().currentMeanEnergy = 0.0;
371 }
372 else {
373 fCache.Get().currentMeanEnergy = weigh
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[nE
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, nAngularP
390 for (i = 0; i < nAngularParameters;
391 theStore.SetCoeff(0, i, theAngular
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(theMan
401 runnin
402 runnin
403 // fill a Legendrestore
404 G4ParticleHPLegendreStore theStore(2
405 theStore.Init(0, e1, nAngularParamet
406 theStore.Init(1, e2, nAngularParamet
407 for (i = 0; i < nAngularParameters;
408 theStore.SetCoeff(0, i, theAngular
409 theStore.SetCoeff(1, i, theAngular
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] / runnin
418 G4double x2 = running[it] / running[nE
419 G4double y1 = theAngular[it - 1].GetLa
420 G4double y2 = theAngular[it].GetLabel(
421 fsEnergy = theInt.Interpolate(theManag
422 G4ParticleHPLegendreStore theStore(2);
423 theStore.Init(0, y1, nAngularParameter
424 theStore.Init(1, y2, nAngularParameter
425 theStore.SetManager(theManager);
426 for (i = 0; i < nAngularParameters; i+
427 theStore.SetCoeff(0, i, theAngular[i
428 theStore.SetCoeff(1, i, theAngular[i
429 }
430 // use it to sample.
431 cosTh = theStore.SampleMax(fsEnergy);
432 }
433 delete[] running;
434
435 // The remaining energy needs to be lowe
436 // *any* case. Otherwise additional pho
437 // produced - therefore the adjustResul
438
439 fCache.Get().remaining_energy -= fsEnerg
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
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(theM
453 theA
454 theA
455 weighted += theInt.GetWeightedBinIntegra
456 theManager.GetScheme(j - 1), theAngula
457 theAngular[j - 1].GetValue(0), theAngu
458 }
459
460 // Cache the mean energy in this distribut
461 if (nEnergies == 1)
462 fCache.Get().currentMeanEnergy = 0.0;
463 else
464 fCache.Get().currentMeanEnergy = weighte
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] < run
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.G
484
485 // For theta, interpolate the compoundFrac
486 G4double cLow = theAngular[itt - 1].GetVal
487 G4double cHigh = theAngular[itt].GetValue(
488 compoundFraction = theInt.Interpolate(theM
489
490 if (compoundFraction > 1.0)
491 compoundFraction = 1.0; // Protection a
492
493 delete[] running;
494
495 // get cosTh
496 G4double incidentEnergy = anEnergy;
497 G4double incidentMass = theProjectile->Get
498 G4double productEnergy = fsEnergy;
499 G4double productMass = result->GetMass();
500 auto targetZ = G4int(fCache.Get().theTarge
501 auto targetA = G4int(fCache.Get().theTarge
502
503 // To correspond to natural composition (-
504 if (targetA == 0) targetA = G4int(fCache.G
505 G4double targetMass = fCache.Get().theTarg
506 auto incidentA = G4int(incidentMass / amu_
507 auto incidentZ = G4int(theProjectile->GetP
508 G4int residualA = targetA + incidentA - A;
509 G4int residualZ = targetZ + incidentZ - Z;
510 G4double residualMass = G4NucleiProperties
511
512 G4ParticleHPKallbachMannSyst theKallbach(
513 compoundFraction, incidentEnergy, incide
514 residualA, residualZ, targetMass, target
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(theM
526 theA
527 theA
528 weighted += theInt.GetWeightedBinIntegra
529 theManager.GetScheme(i - 1), theAngula
530 theAngular[i - 1].GetValue(0), theAngu
531 }
532
533 // Cache the mean energy in this distribut
534 if (nEnergies == 1)
535 fCache.Get().currentMeanEnergy = 0.0;
536 else
537 fCache.Get().currentMeanEnergy = weighte
538
539 if (nEnergies == 1) it = 0;
540 for (i = 1; i < nEnergies; i++) {
541 it = i;
542 if (random < running[i] / running[nEnerg
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 < nAngularParamete
551 theStore.SetX(aCounter, theAngular[0
552 theStore.SetY(aCounter, theAngular[0
553 aCounter++;
554 }
555 G4InterpolationManager aMan;
556 aMan.Init(angularRep - 10, nAngularPar
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, nAngularPar
565 theStore.SetInterpolationManager(aMan)
566 G4InterpolationScheme currentScheme =
567 G4int aCounter = 0;
568 for (G4int j = 1; j < nAngularParamete
569 theStore.SetX(aCounter, theAngular[i
570 theStore.SetY(aCounter, theInt.Inter
571
572
573
574
575 ++aCounter;
576 }
577 cosTh = theStore.Sample();
578 }
579 }
580 else {
581 G4double x1 = running[it - 1] / running[
582 G4double x2 = running[it] / running[nEne
583 G4double y1 = theAngular[it - 1].GetLabe
584 G4double y2 = theAngular[it].GetLabel();
585 fsEnergy = theInt.Interpolate(theManager
586 G4ParticleHPVector theBuff1;
587 G4ParticleHPVector theBuff2;
588 G4InterpolationManager aMan;
589 aMan.Init(angularRep - 10, nAngularParam
590
591 G4int j;
592 for (i = 0, j = 1; i < nAngularParameter
593 theBuff1.SetX(i, theAngular[it - 1].Ge
594 theBuff1.SetY(i, theAngular[it - 1].Ge
595 theBuff2.SetX(i, theAngular[it].GetVal
596 theBuff2.SetY(i, theAngular[it].GetVal
597 }
598
599 G4ParticleHPVector theStore;
600 theStore.SetInterpolationManager(aMan);
601 x1 = y1;
602 x2 = y2;
603 G4double x, y;
604 for (i = 0; i < theBuff1.GetVectorLength
605 x = theBuff1.GetX(i); // costh binnin
606 y1 = theBuff1.GetY(i);
607 y2 = theBuff2.GetY(i);
608 y = theInt.Interpolate(theManager.GetS
609 theAngular[it].
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 "G4ParticleHPCon
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.
628 G4double mtot = result->GetTotalMomentum();
629 G4ThreeVector tempVector(mtot * sinth * std:
630 result->SetMomentum(tempVector);
631 return result;
632 }
633
634 void G4ParticleHPContAngularPar::PrepareTableI
635 {
636 // Discrete energies: store own energies in
637 //
638 // The data files sometimes have identical d
639 // which would lead to overwriting the alrea
640 // creating a hole in the lookup table.
641 // No attempt is made here to correct for th
642 // is subtracted from the energy in order to
643
644 for (G4int ie = 0; ie < nDiscreteEnergies; i
645 // check if energy is already present and
646 G4double myE = theAngular[ie].GetLabel();
647 while (theDiscreteEnergiesOwn.find(myE) !=
648 myE -= 1e-6;
649 }
650 theDiscreteEnergiesOwn[myE] = ie;
651 }
652 return;
653 }
654
655 void G4ParticleHPContAngularPar::BuildByInterp
656
657
658
659 {
660 G4int ie, ie1, ie2, ie1Prev, ie2Prev;
661 // Only rebuild the interpolation table if t
662 // For several subsequent samplings of final
663 // interaction the existing table should be
664 if (!fCache.Get().fresh) return;
665
666 // Make copies of angpar1 and angpar2. Since
667 // it can not be excluded that one of them i
668 // potentially the old "this" for creating t
669 // memory corruption if the old is not store
670 const G4ParticleHPContAngularPar copyAngpar1
671
672 nAngularParameters = copyAngpar1.nAngularPar
673 theManager = copyAngpar1.theManager;
674 theEnergy = anEnergy;
675 theMinEner = DBL_MAX; // min and max will b
676 theMaxEner = -DBL_MAX;
677
678 // The two discrete sets must be merged. A v
679 // be copied to the array in the end. Since
680 // contains pointers, can't simply assign el
681 // needs to call the explicit Set() method i
682
683 // First, average probabilities for those li
684 const std::map<G4double, G4int> discEnerOwn1
685 const std::map<G4double, G4int> discEnerOwn2
686 std::map<G4double, G4int>::const_iterator it
687 std::map<G4double, G4int>::const_iterator it
688 std::vector<G4ParticleHPList*> vAngular(disc
689 G4double discEner1;
690 for (itedeo1 = discEnerOwn1.cbegin(); itedeo
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.theAng
708 G4double val1, val2;
709 for (G4int ip = 0; ip < nAngularParameters
710 val1 = copyAngpar1.theAngular[ie1].GetVa
711 if (ie2 != -1) {
712 val2 = copyAngpar2.theAngular[ie2].Get
713 }
714 else {
715 val2 = 0.;
716 }
717 G4double value = theInt.Interpolate(aSch
718 copy
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_iterat
725 G4double discEner2;
726 for (itedeo2 = discEnerOwn2.cbegin(); itedeo
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 != vAn
746 if (discEner2 > (*itv)->GetLabel()) {
747 itv = vAngular.insert(itv, new G4Par
748 (*itv)->SetLabel(copyAngpar2.theAngu
749 isInserted = true;
750 break;
751 }
752 }
753 if (!isInserted) {
754 ie = (G4int)vAngular.size();
755 vAngular.push_back(new G4ParticleHPLis
756 vAngular[ie]->SetLabel(copyAngpar2.the
757 isInserted = true;
758 }
759
760 G4double val1, val2;
761 for (G4int ip = 0; ip < nAngularParamete
762 val1 = 0;
763 val2 = copyAngpar2.theAngular[ie2].Get
764 G4double value = theInt.Interpolate(aS
765 co
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[nDiscret
777 }
778 theDiscreteEnergiesOwn.clear();
779 theDiscreteEnergies.clear();
780 for (ie = 0; ie < nDiscreteEnergies; ++ie) {
781 theAngular[ie].SetLabel(vAngular[ie]->GetL
782 for (G4int ip = 0; ip < nAngularParameters
783 theAngular[ie].SetValue(ip, vAngular[ie]
784 }
785 theDiscreteEnergiesOwn[theAngular[ie].GetL
786 theDiscreteEnergies.insert(theAngular[ie].
787 }
788
789 // The continuous energies need to be made f
790 // ones. Therefore the re-assignemnt of theA
791 // after the continuous energy set is also f
792 // total number of nEnergies is known and th
793
794 // Get minimum and maximum energy interpolat
795 // Don't use theMinEner or theMaxEner here,
796 // need the interpolated range from the orig
797 G4double interMinEner = copyAngpar1.GetMinEn
798 + (theEnergy - copyA
799 * (copyAngpar2.G
800 / (copyAngpar2.G
801 G4double interMaxEner = copyAngpar1.GetMaxEn
802 + (theEnergy - copyA
803 * (copyAngpar2.G
804 / (copyAngpar2.G
805
806 // Loop to energies of new set
807 theEnergiesTransformed.clear();
808
809 G4int nEnergies1 = copyAngpar1.GetNEnergies(
810 G4int nDiscreteEnergies1 = copyAngpar1.GetND
811 G4double minEner1 = copyAngpar1.GetMinEner()
812 G4double maxEner1 = copyAngpar1.GetMaxEner()
813 G4int nEnergies2 = copyAngpar2.GetNEnergies(
814 G4int nDiscreteEnergies2 = copyAngpar2.GetND
815 G4double minEner2 = copyAngpar2.GetMinEner()
816 G4double maxEner2 = copyAngpar2.GetMaxEner()
817
818 // First build the list of transformed energ
819 // to the new min max by assuming that the m
820 // each set would be scalable to the new, in
821 // max range
822
823 G4double e1(0.);
824 G4double eTNorm1(0.);
825 for (ie1 = nDiscreteEnergies1; ie1 < nEnergi
826 e1 = copyAngpar1.theAngular[ie1].GetLabel(
827 eTNorm1 = (e1 - minEner1);
828 if (maxEner1 != minEner1) eTNorm1 /= (maxE
829 if (eTNorm1 >= 0 && eTNorm1 <= 1) theEnerg
830 }
831
832 G4double e2(0.);
833 G4double eTNorm2(0.);
834 for (ie2 = nDiscreteEnergies2; ie2 < nEnergi
835 e2 = copyAngpar2.theAngular[ie2].GetLabel(
836 eTNorm2 = (e2 - minEner2);
837 if (maxEner2 != minEner2) eTNorm2 /= (maxE
838 if (eTNorm2 >= 0 && eTNorm2 <= 1) theEnerg
839 }
840
841 // Now the list of energies is complete
842 nEnergies = nDiscreteEnergies + (G4int)theEn
843
844 // Create final array of angular parameters
845 const std::size_t esize = nEnergies > 0 ? nE
846 auto theNewAngular = new G4ParticleHPList[es
847
848 // Copy discrete energies and interpolated p
849
850 if (theAngular != nullptr) {
851 for (ie = 0; ie < nDiscreteEnergies; ++ie)
852 theNewAngular[ie].SetLabel(theAngular[ie
853 for (G4int ip = 0; ip < nAngularParamete
854 theNewAngular[ie].SetValue(ip, theAngu
855 }
856 }
857 delete[] theAngular;
858 }
859 theAngular = theNewAngular;
860
861 // Interpolate the continuous energies for n
862 auto iteet = theEnergiesTransformed.begin();
863
864 G4double e1Interp(0.);
865 G4double e2Interp(0.);
866 for (ie = nDiscreteEnergies; ie < nEnergies;
867 G4double eT = (*iteet);
868
869 //--- Use eT1 = eT: Get energy and paramet
870 e1Interp = (maxEner1 - minEner1) * eT + mi
871 //----- Get parameter value corresponding
872 for (ie1 = nDiscreteEnergies1; ie1 < nEner
873 if ((copyAngpar1.theAngular[ie1].GetLabe
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 paramet
883 e2Interp = (maxEner2 - minEner2) * eT + mi
884 //----- Get parameter value corresponding
885 for (ie2 = nDiscreteEnergies2; ie2 < nEner
886 if ((copyAngpar2.theAngular[ie2].GetLabe
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 tran
896 G4double eN = (interMaxEner - 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
910 val1 = theInt.Interpolate2(
911 theManager.GetScheme(ie), e1Int
912 copyAngpar1.theAngular[ie1].Get
913 copyAngpar1.theAngular[ie1].Get
914 * (maxEner1 - minEner1);
915 val2 = theInt.Interpolate2(
916 theManager.GetScheme(ie), e2Int
917 copyAngpar2.theAngular[ie2].Get
918 copyAngpar2.theAngular[ie2].Get
919 * (maxEner2 - minEner2);
920
921 value = theInt.Interpolate(aScheme, anEn
922 val1, val2);
923 if (interMaxEner != interMinEner) {
924 value /= (interMaxEner - interMinEner)
925 }
926 else if (value != 0) {
927 throw G4HadronicException(__FILE__, __
928 "G4ParticleH
929 "interMaxEne
930 }
931 theAngular[ie].SetValue(ip, value);
932 }
933 } // end loop on nDiscreteEnergies
934
935 for (itv = vAngular.cbegin(); itv != vAngula
936 delete (*itv);
937 }
938
939 void G4ParticleHPContAngularPar::Dump() const
940 {
941 G4cout << theEnergy << " " << nEnergies << "
942 << G4endl;
943
944 for (G4int ii = 0; ii < nEnergies; ++ii)
945 theAngular[ii].Dump();
946 }
947