| 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-April-06 Enable IC electron emissions T.
31 // 26-January-07 Add G4NEUTRONHP_USE_ONLY_PHOT
32 // 081024 G4NucleiPropertiesTable:: to G4Nucle
33 // 101203 Bugzilla/Geant4 Problem 1155 Lack of
34 // 110430 Temporary solution in the case of be
35 //
36 // P. Arce, June-2014 Conversion neutron_hp to
37 // V. Ivanchenko July-2023 converted back
38 //
39 #include "G4NeutronHPCaptureFS.hh"
40
41 #include "G4Fragment.hh"
42 #include "G4Gamma.hh"
43 #include "G4IonTable.hh"
44 #include "G4Nucleus.hh"
45 #include "G4ParticleHPDataUsed.hh"
46 #include "G4ParticleHPManager.hh"
47 #include "G4PhotonEvaporation.hh"
48 #include "G4PhysicalConstants.hh"
49 #include "G4PhysicsModelCatalog.hh"
50 #include "G4ReactionProduct.hh"
51 #include "G4RandomDirection.hh"
52 #include "G4SystemOfUnits.hh"
53 #include <sstream>
54
55 G4NeutronHPCaptureFS::G4NeutronHPCaptureFS()
56 {
57 secID = G4PhysicsModelCatalog::GetModelID("m
58 hasXsec = false;
59 hasExactMF6 = false;
60 targetMass = 0;
61 }
62
63 G4HadFinalState*
64 G4NeutronHPCaptureFS::ApplyYourself(const G4Ha
65 {
66 if (theResult.Get() == nullptr) theResult.Pu
67 theResult.Get()->Clear();
68
69 G4int i;
70
71 // prepare neutron
72 G4double eKinetic = theTrack.GetKineticEnerg
73 const G4HadProjectile* incidentParticle = &t
74 G4ReactionProduct theNeutron(theTrack.GetDef
75 theNeutron.SetMomentum(incidentParticle->Get
76 theNeutron.SetKineticEnergy(eKinetic);
77
78 // Prepare target
79 G4ReactionProduct theTarget;
80 G4Nucleus aNucleus;
81 if (targetMass < 500 * MeV)
82 targetMass = G4NucleiProperties::GetNuclea
83 / CLHEP::neutron_mass_c2;
84 G4ThreeVector neutronVelocity = theNeutron.G
85 G4double temperature = theTrack.GetMaterial(
86 theTarget = aNucleus.GetBiasedThermalNucleus
87 theTarget.SetDefinitionAndUpdateE(ionTable->
88
89 // Put neutron in nucleus rest system
90 theNeutron.Lorentz(theNeutron, theTarget);
91 eKinetic = theNeutron.GetKineticEnergy();
92
93 // Sample the photons
94 G4ReactionProductVector* thePhotons = nullpt
95 if (HasFSData() && !fManager->GetUseOnlyPhot
96 // NDL has final state data
97 if (hasExactMF6) {
98 theMF6FinalState.SetTarget(theTarget);
99 theMF6FinalState.SetProjectileRP(theNeut
100 thePhotons = theMF6FinalState.Sample(eKi
101 }
102 else {
103 thePhotons = theFinalStatePhotons.GetPho
104 }
105 if (thePhotons == nullptr) {
106 throw G4HadronicException(__FILE__, __LI
107 "Final state d
108 }
109 }
110 else {
111 // NDL does not have final state data or f
112 G4ThreeVector aCMSMomentum = theNeutron.Ge
113 G4LorentzVector p4(aCMSMomentum, theTarget
114 G4Fragment nucleus(theBaseA + 1, theBaseZ,
115 G4PhotonEvaporation photonEvaporation;
116 // T. K. add
117 photonEvaporation.SetICM(true);
118 G4FragmentVector* products = photonEvapora
119 thePhotons = new G4ReactionProductVector;
120 for (auto it = products->cbegin(); it != p
121 auto theOne = new G4ReactionProduct;
122 // T. K. add
123 if ((*it)->GetParticleDefinition() != nu
124 theOne->SetDefinition((*it)->GetPartic
125 else
126 theOne->SetDefinition(G4Gamma::Gamma()
127
128 // T. K. comment out below line
129 if ((*it)->GetMomentum().mag() > 10 * CL
130 theOne->SetDefinition(ionTable->GetIon
131
132 if ((*it)->GetExcitationEnergy() > 1.0e-
133 G4double ex = (*it)->GetExcitationEner
134 auto aPhoton = new G4ReactionProduct;
135 aPhoton->SetDefinition(G4Gamma::Gamma(
136 aPhoton->SetMomentum((*it)->GetMomentu
137 // aPhoton->SetTotalEnergy( ex ); //wi
138 thePhotons->push_back(aPhoton);
139 }
140
141 theOne->SetMomentum((*it)->GetMomentum()
142 * ((*it)->GetMomentu
143 / (*it)->GetMomentum
144 thePhotons->push_back(theOne);
145 delete *it;
146 }
147 delete products;
148 }
149
150 // Add them to the final state
151 G4int nPhotons = (G4int)thePhotons->size();
152
153 if (!fManager->GetDoNotAdjustFinalState()) {
154 // Make at least one photon
155 // 101203 TK
156 if (nPhotons == 0) {
157 auto theOne = new G4ReactionProduct;
158 theOne->SetDefinition(G4Gamma::Gamma());
159 G4ThreeVector direction = G4RandomDirect
160 theOne->SetMomentum(direction);
161 thePhotons->push_back(theOne);
162 ++nPhotons; // 0 -> 1
163 }
164 // One photon case: energy set to Q-value
165 // 101203 TK
166 if (nPhotons == 1 &&
167 (*thePhotons)[0]->GetDefinition()->GetBaryon
168 G4ThreeVector direction = (*thePhotons)[
169
170 G4double Q = ionTable->GetIonMass(theBas
171 + CLHEP::neutron_mass_c2
172 - ionTable->GetIonMass(theBaseZ, theBaseA +
173 (*thePhotons)[0]->SetMomentum(Q * direct
174 }
175 }
176
177 // back to lab system
178 for (i = 0; i < nPhotons; i++) {
179 (*thePhotons)[i]->Lorentz(*((*thePhotons)[
180 }
181
182 // Recoil, if only one gamma
183 // if (1==nPhotons)
184 if (nPhotons == 1 && thePhotons->operator[](
185 auto theOne = new G4DynamicParticle;
186 G4ParticleDefinition* aRecoil = ionTable->
187 theOne->SetDefinition(aRecoil);
188 // Now energy;
189 // Can be done slightly better @
190 G4ThreeVector aMomentum = theTrack.Get4Mom
191 - thePhotons->op
192
193 theOne->SetMomentum(aMomentum);
194 theResult.Get()->AddSecondary(theOne, secI
195 }
196
197 // Now fill in the gammas.
198 for (i = 0; i < nPhotons; ++i) {
199 // back to lab system
200 auto theOne = new G4DynamicParticle;
201 theOne->SetDefinition(thePhotons->operator
202 theOne->SetMomentum(thePhotons->operator[]
203 theResult.Get()->AddSecondary(theOne, secI
204 delete thePhotons->operator[](i);
205 }
206 delete thePhotons;
207
208 // 101203TK
209 G4bool residual = false;
210 G4ParticleDefinition* aRecoil = ionTable->Ge
211 for (std::size_t j = 0; j != theResult.Get()
212 if (theResult.Get()->GetSecondary(j)->GetP
213 residual = true;
214 }
215
216 if (!residual) {
217 G4int nNonZero = 0;
218 G4LorentzVector p_photons(0, 0, 0, 0);
219 for (std::size_t j = 0; j != theResult.Get
220 p_photons += theResult.Get()->GetSeconda
221 // To many 0 momentum photons -> Check P
222 if (theResult.Get()->GetSecondary(j)->Ge
223 }
224
225 // Can we include kinetic energy here?
226 G4double deltaE = (theTrack.Get4Momentum()
227 - (p_photons.e() + aReco
228
229 // Add photons
230 if (nPhotons - nNonZero > 0) {
231 // G4cout << "TKDB G4NeutronHPCaptureFS:
232 // nPhotons - nNonZero << " photons" <<
233 std::vector<G4double> vRand;
234 vRand.push_back(0.0);
235 for (G4int j = 0; j != nPhotons - nNonZe
236 vRand.push_back(G4UniformRand());
237 }
238 vRand.push_back(1.0);
239 std::sort(vRand.begin(), vRand.end());
240
241 std::vector<G4double> vEPhoton;
242 for (G4int j = 0; j < (G4int)vRand.size(
243 vEPhoton.push_back(deltaE * (vRand[j +
244 }
245 std::sort(vEPhoton.begin(), vEPhoton.end
246
247 for (G4int j = 0; j < nPhotons - nNonZer
248 // Isotopic in LAB OK?
249 // Bug # 1745 DHW G4double theta = pi
250 G4ThreeVector tempVector = G4RandomDir
251
252 p_photons += G4LorentzVector(tempVecto
253 auto theOne = new G4DynamicParticle;
254 theOne->SetDefinition(G4Gamma::Gamma()
255 theOne->SetMomentum(tempVector);
256 theResult.Get()->AddSecondary(theOne,
257 }
258
259 // Add last photon
260 auto theOne = new G4DynamicParticle;
261 theOne->SetDefinition(G4Gamma::Gamma());
262 // For better momentum conservati
263 G4ThreeVector lastPhoton = -p_photons.ve
264 p_photons += G4LorentzVector(lastPhoton,
265 theOne->SetMomentum(lastPhoton);
266 theResult.Get()->AddSecondary(theOne, se
267 }
268
269 // Add residual
270 auto theOne = new G4DynamicParticle;
271 G4ThreeVector aMomentum =
272 theTrack.Get4Momentum().vect() + theTarg
273 theOne->SetDefinition(aRecoil);
274 theOne->SetMomentum(aMomentum);
275 theResult.Get()->AddSecondary(theOne, secI
276 }
277 // 101203TK END
278
279 // clean up the primary neutron
280 theResult.Get()->SetStatusChange(stopAndKill
281 return theResult.Get();
282 }
283
284 void G4NeutronHPCaptureFS::Init(G4double AA, G
285 const G4String
286 G4ParticleDefi
287 {
288 G4int Z = G4lrint(ZZ);
289 G4int A = G4lrint(AA);
290 // TK110430 BEGIN
291 std::stringstream ss;
292 ss << Z;
293 G4String sZ;
294 ss >> sZ;
295 ss.clear();
296 ss << A;
297 G4String sA;
298 ss >> sA;
299
300 ss.clear();
301 G4String sM;
302 if (M > 0) {
303 ss << "m";
304 ss << M;
305 ss >> sM;
306 ss.clear();
307 }
308
309 G4String element_name = theNames.GetName(Z -
310 G4String filenameMF6 = dirName + "/FSMF6/" +
311
312 std::istringstream theData(std::ios::in);
313 G4ParticleHPManager::GetInstance()->GetDataS
314
315 // TK110430 Only use MF6MT102 which has exac
316 // Even _nat_ do not select and there is no
317 if (theData.good()) {
318 hasExactMF6 = true;
319 theMF6FinalState.Init(theData);
320 return;
321 }
322 // TK110430 END
323
324 const G4String& tString = "/FS";
325 G4bool dbool;
326 const G4ParticleHPDataUsed& aFile =
327 theNames.GetName(A, Z, M, dirName, tString
328
329 const G4String& filename = aFile.GetName();
330 SetAZMs(A, Z, M, aFile);
331 if (!dbool || (Z <= 2 && (theBaseZ != Z || t
332 hasAnyData = false;
333 hasFSData = false;
334 hasXsec = false;
335 return;
336 }
337 theData.clear();
338 G4ParticleHPManager::GetInstance()->GetDataS
339 hasFSData = theFinalStatePhotons.InitMean(th
340 if (hasFSData) {
341 targetMass = theFinalStatePhotons.GetTarge
342 theFinalStatePhotons.InitAngular(theData);
343 theFinalStatePhotons.InitEnergies(theData)
344 }
345 }
346