| 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 // 080801 Protect div0 error, when theCompundF
31 //
32 // P. Arce, June-2014 Conversion neutron_hp to
33 //
34 // June-2019 - E. Mendoza --> perform some cor
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(
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 thr
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,
72
73 return result;
74 }
75
76 G4double G4ParticleHPKallbachMannSyst::Kallbac
77 {
78 // Kallbach-Mann systematics without normali
79 G4double result;
80 G4double theX = A(anEnergy) * cosTh;
81 result =
82 0.5 * (G4Exp(theX) * (1 + theCompoundFract
83 return result;
84 }
85
86 G4double G4ParticleHPKallbachMannSyst::GetKall
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) < delt
94 theCompoundFraction = 1.0 - delta;
95 }
96 result = 0.5 * (1. / A(anEnergy)) * G4Log((1
97 return result;
98 }
99
100 G4double G4ParticleHPKallbachMannSyst::A(G4dou
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 * Me
106
107 G4double epsa = anEnergy * theTargetMass / (
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
113 G4double Et1 = 130 * MeV;
114 G4double R1 = std::min(ea, Et1);
115 // theProductEnergy is still in CMS!!!
116 G4double epsb = theProductEnergy * (theProdu
117 G4int AB = theResidualA;
118 G4int ZB = theResidualZ;
119 G4double eb = epsb + SeparationEnergy(Ac, Nc
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 =
128 Ma = 1;
129 } // neutron,proton,deuteron
130 else if (theProjectileA == 4 && theProjectil
131 Ma = 0;
132 } // alpha
133 else if (theProjectileA == 3 && (theProjecti
134 Ma = 0.5;
135 } // tritum,He3 : set intermediate value
136 else {
137 throw G4HadronicException(__FILE__, __LINE
138 "Severe error in
139 }
140 if (theProductA == 1 && theProductZ == 0) {
141 mb = 1. / 2.;
142 } // neutron
143 else if (theProductA == 4 && theProductZ ==
144 mb = 2;
145 } // alpha
146 else {
147 mb = 1;
148 }
149
150 result = C1 * X1 + C2 * G4Pow::GetInstance()
151 + C3 * Ma * mb * G4Pow::GetInstance
152 return result;
153 }
154
155 G4double G4ParticleHPKallbachMannSyst::Separat
156
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) /
163 result +=
164 -18.56 * (G4Pow::GetInstance()->A23(G4doub
165 result += 33.22
166 * ((Nc - Zc) * (Nc - Zc) / G4Pow::
167 - (NA - ZA) * (NA - ZA) / G4Pow
168 result += -0.717
169 * (Zc * Zc / G4Pow::GetInstance()-
170 - ZA * ZA / G4Pow::GetInstance(
171 result += 1.211 * (Zc * Zc / (G4double)Ac -
172 G4double totalBinding(0);
173 if (Zbinding == 0 && Abinding == 1) totalBin
174 if (Zbinding == 1 && Abinding == 1) totalBin
175 if (Zbinding == 1 && Abinding == 2) totalBin
176 if (Zbinding == 1 && Abinding == 3) totalBin
177 if (Zbinding == 2 && Abinding == 3) totalBin
178 if (Zbinding == 2 && Abinding == 4) totalBin
179 result += -totalBinding;
180 result *= MeV;
181 return result;
182 }
183