| 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 // 12-April-06 Enable IC electron emissions T. Koi
31 // 26-January-07 Add G4NEUTRONHP_USE_ONLY_PHOTONEVAPORATION flag
32 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
33 // 101203 Bugzilla/Geant4 Problem 1155 Lack of residual in some case
34 // 110430 Temporary solution in the case of being MF6 final state in Capture reaction (MT102)
35 //
36 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
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("model_NeutronHPCapture");
58 hasXsec = false;
59 hasExactMF6 = false;
60 targetMass = 0;
61 }
62
63 G4HadFinalState*
64 G4NeutronHPCaptureFS::ApplyYourself(const G4HadProjectile& theTrack)
65 {
66 if (theResult.Get() == nullptr) theResult.Put(new G4HadFinalState);
67 theResult.Get()->Clear();
68
69 G4int i;
70
71 // prepare neutron
72 G4double eKinetic = theTrack.GetKineticEnergy();
73 const G4HadProjectile* incidentParticle = &theTrack;
74 G4ReactionProduct theNeutron(theTrack.GetDefinition());
75 theNeutron.SetMomentum(incidentParticle->Get4Momentum().vect());
76 theNeutron.SetKineticEnergy(eKinetic);
77
78 // Prepare target
79 G4ReactionProduct theTarget;
80 G4Nucleus aNucleus;
81 if (targetMass < 500 * MeV)
82 targetMass = G4NucleiProperties::GetNuclearMass(theBaseA, theBaseZ)
83 / CLHEP::neutron_mass_c2;
84 G4ThreeVector neutronVelocity = theNeutron.GetMomentum()/ CLHEP::neutron_mass_c2;
85 G4double temperature = theTrack.GetMaterial()->GetTemperature();
86 theTarget = aNucleus.GetBiasedThermalNucleus(targetMass, neutronVelocity, temperature);
87 theTarget.SetDefinitionAndUpdateE(ionTable->GetIon(theBaseZ, theBaseA, 0.0));
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 = nullptr;
95 if (HasFSData() && !fManager->GetUseOnlyPhotoEvaporation()) {
96 // NDL has final state data
97 if (hasExactMF6) {
98 theMF6FinalState.SetTarget(theTarget);
99 theMF6FinalState.SetProjectileRP(theNeutron);
100 thePhotons = theMF6FinalState.Sample(eKinetic);
101 }
102 else {
103 thePhotons = theFinalStatePhotons.GetPhotons(eKinetic);
104 }
105 if (thePhotons == nullptr) {
106 throw G4HadronicException(__FILE__, __LINE__,
107 "Final state data for photon is not properly allocated");
108 }
109 }
110 else {
111 // NDL does not have final state data or forced to use PhotoEvaporation model
112 G4ThreeVector aCMSMomentum = theNeutron.GetMomentum() + theTarget.GetMomentum();
113 G4LorentzVector p4(aCMSMomentum, theTarget.GetTotalEnergy() + theNeutron.GetTotalEnergy());
114 G4Fragment nucleus(theBaseA + 1, theBaseZ, p4);
115 G4PhotonEvaporation photonEvaporation;
116 // T. K. add
117 photonEvaporation.SetICM(true);
118 G4FragmentVector* products = photonEvaporation.BreakItUp(nucleus);
119 thePhotons = new G4ReactionProductVector;
120 for (auto it = products->cbegin(); it != products->cend(); ++it) {
121 auto theOne = new G4ReactionProduct;
122 // T. K. add
123 if ((*it)->GetParticleDefinition() != nullptr)
124 theOne->SetDefinition((*it)->GetParticleDefinition());
125 else
126 theOne->SetDefinition(G4Gamma::Gamma()); // this definiion will be over writen
127
128 // T. K. comment out below line
129 if ((*it)->GetMomentum().mag() > 10 * CLHEP::MeV)
130 theOne->SetDefinition(ionTable->GetIon(theBaseZ, theBaseA + 1, 0));
131
132 if ((*it)->GetExcitationEnergy() > 1.0e-2 * eV) {
133 G4double ex = (*it)->GetExcitationEnergy();
134 auto aPhoton = new G4ReactionProduct;
135 aPhoton->SetDefinition(G4Gamma::Gamma());
136 aPhoton->SetMomentum((*it)->GetMomentum().vect().unit() * ex);
137 // aPhoton->SetTotalEnergy( ex ); //will be calculated from momentum
138 thePhotons->push_back(aPhoton);
139 }
140
141 theOne->SetMomentum((*it)->GetMomentum().vect()
142 * ((*it)->GetMomentum().t() - (*it)->GetExcitationEnergy())
143 / (*it)->GetMomentum().t());
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 = G4RandomDirection();
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()->GetBaryonNumber() == 0) {
168 G4ThreeVector direction = (*thePhotons)[0]->GetMomentum().unit();
169
170 G4double Q = ionTable->GetIonMass(theBaseZ, theBaseA, 0)
171 + CLHEP::neutron_mass_c2
172 - ionTable->GetIonMass(theBaseZ, theBaseA + 1, 0);
173 (*thePhotons)[0]->SetMomentum(Q * direction);
174 }
175 }
176
177 // back to lab system
178 for (i = 0; i < nPhotons; i++) {
179 (*thePhotons)[i]->Lorentz(*((*thePhotons)[i]), -1*theTarget);
180 }
181
182 // Recoil, if only one gamma
183 // if (1==nPhotons)
184 if (nPhotons == 1 && thePhotons->operator[](0)->GetDefinition()->GetBaryonNumber() == 0) {
185 auto theOne = new G4DynamicParticle;
186 G4ParticleDefinition* aRecoil = ionTable->GetIon(theBaseZ, theBaseA + 1, 0);
187 theOne->SetDefinition(aRecoil);
188 // Now energy;
189 // Can be done slightly better @
190 G4ThreeVector aMomentum = theTrack.Get4Momentum().vect() + theTarget.GetMomentum()
191 - thePhotons->operator[](0)->GetMomentum();
192
193 theOne->SetMomentum(aMomentum);
194 theResult.Get()->AddSecondary(theOne, secID);
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[](i)->GetDefinition());
202 theOne->SetMomentum(thePhotons->operator[](i)->GetMomentum());
203 theResult.Get()->AddSecondary(theOne, secID);
204 delete thePhotons->operator[](i);
205 }
206 delete thePhotons;
207
208 // 101203TK
209 G4bool residual = false;
210 G4ParticleDefinition* aRecoil = ionTable->GetIon(theBaseZ, theBaseA + 1, 0);
211 for (std::size_t j = 0; j != theResult.Get()->GetNumberOfSecondaries(); j++) {
212 if (theResult.Get()->GetSecondary(j)->GetParticle()->GetDefinition() == aRecoil)
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()->GetNumberOfSecondaries(); ++j) {
220 p_photons += theResult.Get()->GetSecondary(j)->GetParticle()->Get4Momentum();
221 // To many 0 momentum photons -> Check PhotonDist
222 if (theResult.Get()->GetSecondary(j)->GetParticle()->Get4Momentum().e() > 0) nNonZero++;
223 }
224
225 // Can we include kinetic energy here?
226 G4double deltaE = (theTrack.Get4Momentum().e() + theTarget.GetTotalEnergy())
227 - (p_photons.e() + aRecoil->GetPDGMass());
228
229 // Add photons
230 if (nPhotons - nNonZero > 0) {
231 // G4cout << "TKDB G4NeutronHPCaptureFS::ApplyYourself we will create additional " <<
232 // nPhotons - nNonZero << " photons" << G4endl;
233 std::vector<G4double> vRand;
234 vRand.push_back(0.0);
235 for (G4int j = 0; j != nPhotons - nNonZero - 1; j++) {
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() - 1; j++) {
243 vEPhoton.push_back(deltaE * (vRand[j + 1] - vRand[j]));
244 }
245 std::sort(vEPhoton.begin(), vEPhoton.end());
246
247 for (G4int j = 0; j < nPhotons - nNonZero - 1; j++) {
248 // Isotopic in LAB OK?
249 // Bug # 1745 DHW G4double theta = pi*G4UniformRand();
250 G4ThreeVector tempVector = G4RandomDirection()*vEPhoton[j];
251
252 p_photons += G4LorentzVector(tempVector, tempVector.mag());
253 auto theOne = new G4DynamicParticle;
254 theOne->SetDefinition(G4Gamma::Gamma());
255 theOne->SetMomentum(tempVector);
256 theResult.Get()->AddSecondary(theOne, secID);
257 }
258
259 // Add last photon
260 auto theOne = new G4DynamicParticle;
261 theOne->SetDefinition(G4Gamma::Gamma());
262 // For better momentum conservation
263 G4ThreeVector lastPhoton = -p_photons.vect().unit() * vEPhoton.back();
264 p_photons += G4LorentzVector(lastPhoton, lastPhoton.mag());
265 theOne->SetMomentum(lastPhoton);
266 theResult.Get()->AddSecondary(theOne, secID);
267 }
268
269 // Add residual
270 auto theOne = new G4DynamicParticle;
271 G4ThreeVector aMomentum =
272 theTrack.Get4Momentum().vect() + theTarget.GetMomentum() - p_photons.vect();
273 theOne->SetDefinition(aRecoil);
274 theOne->SetMomentum(aMomentum);
275 theResult.Get()->AddSecondary(theOne, secID);
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, G4double ZZ, G4int M,
285 const G4String& dirName, const G4String&,
286 G4ParticleDefinition*)
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 - 1);
310 G4String filenameMF6 = dirName + "/FSMF6/" + sZ + "_" + sA + sM + "_" + element_name;
311
312 std::istringstream theData(std::ios::in);
313 G4ParticleHPManager::GetInstance()->GetDataStream(filenameMF6, theData);
314
315 // TK110430 Only use MF6MT102 which has exactly same A and Z
316 // Even _nat_ do not select and there is no _nat_ case in ENDF-VII.0
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, dbool);
328
329 const G4String& filename = aFile.GetName();
330 SetAZMs(A, Z, M, aFile);
331 if (!dbool || (Z <= 2 && (theBaseZ != Z || theBaseA != A))) {
332 hasAnyData = false;
333 hasFSData = false;
334 hasXsec = false;
335 return;
336 }
337 theData.clear();
338 G4ParticleHPManager::GetInstance()->GetDataStream(filename, theData);
339 hasFSData = theFinalStatePhotons.InitMean(theData);
340 if (hasFSData) {
341 targetMass = theFinalStatePhotons.GetTargetMass();
342 theFinalStatePhotons.InitAngular(theData);
343 theFinalStatePhotons.InitEnergies(theData);
344 }
345 }
346