| 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 // 12-Apr-06 fix in delayed neutron and photon
31 // 07-Sep-11 M. Kelsey -- Follow change to G4H
32 // P. Arce, June-2014 Conversion neutron_hp to
33 //
34
35 #include "G4ParticleHPFissionFS.hh"
36
37 #include "G4DynamicParticleVector.hh"
38 #include "G4Exp.hh"
39 #include "G4IonTable.hh"
40 #include "G4Nucleus.hh"
41 #include "G4ParticleHPFissionERelease.hh"
42 #include "G4PhysicalConstants.hh"
43 #include "G4PhysicsModelCatalog.hh"
44
45 G4ParticleHPFissionFS::G4ParticleHPFissionFS()
46 {
47 secID = G4PhysicsModelCatalog::GetModelID("m
48 hasXsec = false;
49 produceFissionFragments = false;
50 }
51
52 void G4ParticleHPFissionFS::Init(G4double A, G
53 const G4Strin
54 {
55 theFS.Init(A, Z, M, dirName, aFSType, projec
56 theFC.Init(A, Z, M, dirName, aFSType, projec
57 theSC.Init(A, Z, M, dirName, aFSType, projec
58 theTC.Init(A, Z, M, dirName, aFSType, projec
59 theLC.Init(A, Z, M, dirName, aFSType, projec
60
61 theFF.Init(A, Z, M, dirName, aFSType, projec
62 if (G4ParticleHPManager::GetInstance()->GetP
63 G4cout << "Fission fragment production is
64 << "Z = " << (G4int)Z << ", A = " <
65 G4cout << "As currently modeled this optio
66 "fission fragments."
67 << G4endl;
68 produceFissionFragments = true;
69 }
70 }
71
72 G4HadFinalState* G4ParticleHPFissionFS::ApplyY
73 {
74 // Because it may change by UI command
75 produceFissionFragments = G4ParticleHPManage
76
77 // prepare neutron
78 if (theResult.Get() == nullptr) theResult.Pu
79 theResult.Get()->Clear();
80 G4double eKinetic = theTrack.GetKineticEnerg
81 const G4HadProjectile* incidentParticle = &t
82 G4ReactionProduct theNeutron(
83 const_cast<G4ParticleDefinition*>(incident
84 theNeutron.SetMomentum(incidentParticle->Get
85 theNeutron.SetKineticEnergy(eKinetic);
86
87 // prepare target
88 G4Nucleus aNucleus;
89 G4ReactionProduct theTarget;
90 G4double targetMass = theFS.GetMass();
91 G4ThreeVector neuVelo =
92 (1. / incidentParticle->GetDefinition()->G
93 theTarget =
94 aNucleus.GetBiasedThermalNucleus(targetMas
95 theTarget.SetDefinition(
96 G4IonTable::GetIonTable()->GetIon(G4int(th
97
98 // set neutron and target in the FS classes
99 theFS.SetNeutronRP(theNeutron);
100 theFS.SetTarget(theTarget);
101 theFC.SetNeutronRP(theNeutron);
102 theFC.SetTarget(theTarget);
103 theSC.SetNeutronRP(theNeutron);
104 theSC.SetTarget(theTarget);
105 theTC.SetNeutronRP(theNeutron);
106 theTC.SetTarget(theTarget);
107 theLC.SetNeutronRP(theNeutron);
108 theLC.SetTarget(theTarget);
109 theFF.SetNeutronRP(theNeutron);
110 theFF.SetTarget(theTarget);
111
112 // boost to target rest system and decide on
113 theNeutron.Lorentz(theNeutron, -1 * theTarge
114
115 // dice the photons
116
117 G4DynamicParticleVector* thePhotons;
118 thePhotons = theFS.GetPhotons();
119
120 // select the FS in charge
121
122 eKinetic = theNeutron.GetKineticEnergy();
123 G4double xSec[4];
124 xSec[0] = theFC.GetXsec(eKinetic);
125 xSec[1] = xSec[0] + theSC.GetXsec(eKinetic);
126 xSec[2] = xSec[1] + theTC.GetXsec(eKinetic);
127 xSec[3] = xSec[2] + theLC.GetXsec(eKinetic);
128 G4int it;
129 unsigned int i = 0;
130 G4double random = G4UniformRand();
131 if (xSec[3] == 0) {
132 it = -1;
133 }
134 else {
135 for (i = 0; i < 4; i++) {
136 it = i;
137 if (random < xSec[i] / xSec[3]) break;
138 }
139 }
140
141 // dice neutron multiplicities, energies and
142 // no energy conservation on an event-to-eve
143 // also for mean, we rely on the consistancy
144
145 G4int Prompt = 0, delayed = 0, all = 0;
146 G4DynamicParticleVector* theNeutrons = nullp
147 switch (it) // check logic, and ask, if par
148 // particles @@@
149 {
150 case 0:
151 theFS.SampleNeutronMult(all, Prompt, del
152 if (Prompt == 0 && delayed == 0) Prompt
153 theNeutrons = theFC.ApplyYourself(Prompt
154 // take 'U' into account explicitly (see
155 break;
156 case 1:
157 theFS.SampleNeutronMult(all, Prompt, del
158 if (Prompt == 0 && delayed == 0) Prompt
159 theNeutrons = theSC.ApplyYourself(Prompt
160 break;
161 case 2:
162 theFS.SampleNeutronMult(all, Prompt, del
163 if (Prompt == 0 && delayed == 0) Prompt
164 theNeutrons = theTC.ApplyYourself(Prompt
165 break;
166 case 3:
167 theFS.SampleNeutronMult(all, Prompt, del
168 if (Prompt == 0 && delayed == 0) Prompt
169 theNeutrons = theLC.ApplyYourself(Prompt
170 break;
171 default:
172 break;
173 }
174
175 // dice delayed neutrons and photons, and fa
176 // for Prompt in case channel had no FS data
177
178 if (produceFissionFragments) delayed = 0;
179
180 G4double* theDecayConstants;
181
182 if (theNeutrons != nullptr) {
183 theDecayConstants = new G4double[delayed];
184 for (i = 0; i < theNeutrons->size(); ++i)
185 theResult.Get()->AddSecondary(theNeutron
186 }
187 delete theNeutrons;
188
189 G4DynamicParticleVector* theDelayed = null
190 theDelayed = theFS.ApplyYourself(0, delaye
191 for (i = 0; i < theDelayed->size(); i++) {
192 G4double time = -G4Log(G4UniformRand())
193 time += theTrack.GetGlobalTime();
194 theResult.Get()->AddSecondary(theDelayed
195 theResult.Get()->GetSecondary(theResult.
196 }
197 delete theDelayed;
198 }
199 else {
200 theFS.SampleNeutronMult(all, Prompt, delay
201 theDecayConstants = new G4double[delayed];
202 if (Prompt == 0 && delayed == 0) Prompt =
203 theNeutrons = theFS.ApplyYourself(Prompt,
204 G4int i0;
205 for (i0 = 0; i0 < Prompt; ++i0) {
206 theResult.Get()->AddSecondary(theNeutron
207 }
208
209 for (i0 = Prompt; i0 < Prompt + delayed; +
210 // Protect against the very rare case of
211 G4double time = 0.0;
212 if (theDecayConstants[i0 - Prompt] > 1.0
213 time = -G4Log(G4UniformRand()) / theDe
214 }
215 else {
216 G4ExceptionDescription ed;
217 ed << " theDecayConstants[i0-Prompt]="
218 << " -> cannot sample the time :
219 G4Exception("G4ParticleHPFissionFS::Ap
220 }
221
222 time += theTrack.GetGlobalTime();
223 theResult.Get()->AddSecondary(theNeutron
224 theResult.Get()->GetSecondary(theResult.
225 }
226 delete theNeutrons;
227 }
228 delete[] theDecayConstants;
229
230 std::size_t nPhotons = 0;
231 if (thePhotons != nullptr) {
232 nPhotons = thePhotons->size();
233 for (i = 0; i < nPhotons; ++i) {
234 theResult.Get()->AddSecondary(thePhotons
235 }
236 delete thePhotons;
237 }
238
239 // finally deal with local energy deposition
240
241 G4ParticleHPFissionERelease* theERelease = t
242 G4double eDepByFragments = theERelease->GetF
243 // theResult.SetLocalEnergyDeposit(eDepByFra
244 if (!produceFissionFragments) theResult.Get(
245 // clean up the primary neutron
246 theResult.Get()->SetStatusChange(stopAndKill
247
248 if (produceFissionFragments) {
249 G4int fragA_Z = 0;
250 G4int fragA_A = 0;
251 G4int fragA_M = 0;
252 // System is traget rest!
253 theFF.GetAFissionFragment(eKinetic, fragA_
254 if (0 == fragA_A) { return theResult.Get()
255 G4int fragB_Z = (G4int)theBaseZ - fragA_Z;
256 G4int fragB_A = (G4int)theBaseA - fragA_A
257
258 G4IonTable* pt = G4IonTable::GetIonTable()
259 // Excitation energy is not taken into acc
260 G4ParticleDefinition* pdA = pt->GetIon(fra
261 G4ParticleDefinition* pdB = pt->GetIon(fra
262
263 // Isotropic Distribution
264 G4double phi = twopi * G4UniformRand();
265 // Bug #1745 DHW G4double theta = pi*G4Uni
266 G4double costheta = 2. * G4UniformRand() -
267 G4double theta = std::acos(costheta);
268 G4double sinth = std::sin(theta);
269 G4ThreeVector direction(sinth * std::cos(p
270
271 // Just use ENDF value for this
272 G4double ER = eDepByFragments;
273 G4double ma = pdA->GetPDGMass();
274 G4double mb = pdB->GetPDGMass();
275 G4double EA = ER / (1 + ma / mb);
276 G4double EB = ER - EA;
277 auto dpA = new G4DynamicParticle(pdA, dire
278 auto dpB = new G4DynamicParticle(pdB, -dir
279 theResult.Get()->AddSecondary(dpA, secID);
280 theResult.Get()->AddSecondary(dpB, secID);
281 }
282
283 return theResult.Get();
284 }
285