| 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 // 080801 Protect div0 error, when theCompundFraction is 1 by T. Koi
31 //
32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
33 //
34 // June-2019 - E. Mendoza --> perform some corrections
35
36 #include "G4ParticleHPKallbachMannSyst.hh"
37
38 #include "G4Exp.hh"
39 #include "G4HadronicException.hh"
40 #include "G4Log.hh"
41 #include "G4Pow.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "Randomize.hh"
44
45 G4double G4ParticleHPKallbachMannSyst::Sample(G4double anEnergy)
46 {
47 G4double result = 0.;
48
49 G4double zero = GetKallbachZero(anEnergy);
50 if (zero > 1) zero = 1.;
51 if (zero < -1) zero = -1.;
52 G4double max = Kallbach(zero, anEnergy);
53 G4double upper = Kallbach(1., anEnergy);
54 G4double lower = Kallbach(-1., anEnergy);
55 if (upper > max) max = upper;
56 if (lower > max) max = lower;
57 G4double value, random;
58
59 G4int icounter = 0;
60 G4int icounter_max = 1024;
61 do {
62 icounter++;
63 if (icounter > icounter_max) {
64 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of "
65 << __FILE__ << "." << G4endl;
66 break;
67 }
68 result = 2. * G4UniformRand() - 1;
69 value = Kallbach(result, anEnergy) / max;
70 random = G4UniformRand();
71 } while (random > value); // Loop checking, 11.05.2015, T. Koi
72
73 return result;
74 }
75
76 G4double G4ParticleHPKallbachMannSyst::Kallbach(G4double cosTh, G4double anEnergy)
77 {
78 // Kallbach-Mann systematics without normalization.
79 G4double result;
80 G4double theX = A(anEnergy) * cosTh;
81 result =
82 0.5 * (G4Exp(theX) * (1 + theCompoundFraction) + G4Exp(-theX) * (1 - theCompoundFraction));
83 return result;
84 }
85
86 G4double G4ParticleHPKallbachMannSyst::GetKallbachZero(G4double anEnergy)
87 {
88 G4double result;
89 // delta 2.0e-16 in not good.
90 // delta 4.0e-16 is OK
91 // safety factor of 2
92 G4double delta = 8.0e-16;
93 if (std::abs(theCompoundFraction - 1) < delta) {
94 theCompoundFraction = 1.0 - delta;
95 }
96 result = 0.5 * (1. / A(anEnergy)) * G4Log((1 - theCompoundFraction) / (1 + theCompoundFraction));
97 return result;
98 }
99
100 G4double G4ParticleHPKallbachMannSyst::A(G4double anEnergy)
101 {
102 G4double result;
103 G4double C1 = 0.04 / MeV;
104 G4double C2 = 1.8E-6 / (MeV * MeV * MeV);
105 G4double C3 = 6.7E-7 / (MeV * MeV * MeV * MeV);
106
107 G4double epsa = anEnergy * theTargetMass / (theTargetMass + theIncidentMass);
108 G4int Ac = theTargetA + theProjectileA;
109 G4int Nc = Ac - theTargetZ - theProjectileZ;
110 G4int AA = theTargetA;
111 G4int ZA = theTargetZ;
112 G4double ea = epsa + SeparationEnergy(Ac, Nc, AA, ZA, theProjectileA, theProjectileZ);
113 G4double Et1 = 130 * MeV;
114 G4double R1 = std::min(ea, Et1);
115 // theProductEnergy is still in CMS!!!
116 G4double epsb = theProductEnergy * (theProductMass + theResidualMass) / theResidualMass;
117 G4int AB = theResidualA;
118 G4int ZB = theResidualZ;
119 G4double eb = epsb + SeparationEnergy(Ac, Nc, AB, ZB, theProductA, theProductZ);
120 G4double X1 = R1 * eb / ea;
121 G4double Et3 = 41 * MeV;
122 G4double R3 = std::min(ea, Et3);
123 G4double X3 = R3 * eb / ea;
124
125 G4double Ma = 1;
126 G4double mb = 1;
127 if (theProjectileA == 1 || (theProjectileZ == 1 && theProjectileA == 2)) {
128 Ma = 1;
129 } // neutron,proton,deuteron
130 else if (theProjectileA == 4 && theProjectileZ == 2) {
131 Ma = 0;
132 } // alpha
133 else if (theProjectileA == 3 && (theProjectileZ == 1 || theProjectileZ == 2)) {
134 Ma = 0.5;
135 } // tritum,He3 : set intermediate value
136 else {
137 throw G4HadronicException(__FILE__, __LINE__,
138 "Severe error in the sampling of Kallbach-Mann Systematics");
139 }
140 if (theProductA == 1 && theProductZ == 0) {
141 mb = 1. / 2.;
142 } // neutron
143 else if (theProductA == 4 && theProductZ == 2) {
144 mb = 2;
145 } // alpha
146 else {
147 mb = 1;
148 }
149
150 result = C1 * X1 + C2 * G4Pow::GetInstance()->powN(X1, 3)
151 + C3 * Ma * mb * G4Pow::GetInstance()->powN(X3, 4);
152 return result;
153 }
154
155 G4double G4ParticleHPKallbachMannSyst::SeparationEnergy(G4int Ac, G4int Nc, G4int AA, G4int ZA,
156 G4int Abinding, G4int Zbinding)
157 {
158 G4double result;
159 G4int NA = AA - ZA;
160 G4int Zc = Ac - Nc;
161 result = 15.68 * (Ac - AA);
162 result += -28.07 * ((Nc - Zc) * (Nc - Zc) / (G4double)Ac - (NA - ZA) * (NA - ZA) / (G4double)AA);
163 result +=
164 -18.56 * (G4Pow::GetInstance()->A23(G4double(Ac)) - G4Pow::GetInstance()->A23(G4double(AA)));
165 result += 33.22
166 * ((Nc - Zc) * (Nc - Zc) / G4Pow::GetInstance()->powA(G4double(Ac), 4. / 3.)
167 - (NA - ZA) * (NA - ZA) / G4Pow::GetInstance()->powA(G4double(AA), 4. / 3.));
168 result += -0.717
169 * (Zc * Zc / G4Pow::GetInstance()->A13(G4double(Ac))
170 - ZA * ZA / G4Pow::GetInstance()->A13(G4double(AA)));
171 result += 1.211 * (Zc * Zc / (G4double)Ac - ZA * ZA / (G4double)AA);
172 G4double totalBinding(0);
173 if (Zbinding == 0 && Abinding == 1) totalBinding = 0;
174 if (Zbinding == 1 && Abinding == 1) totalBinding = 0;
175 if (Zbinding == 1 && Abinding == 2) totalBinding = 2.224596;
176 if (Zbinding == 1 && Abinding == 3) totalBinding = 8.481798;
177 if (Zbinding == 2 && Abinding == 3) totalBinding = 7.718043;
178 if (Zbinding == 2 && Abinding == 4) totalBinding = 28.29566;
179 result += -totalBinding;
180 result *= MeV;
181 return result;
182 }
183