| 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 // particle_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron transport model.
29 //
30 // 080801 Give a warning message for irregular mass value in data file by T. Koi
31 // Introduce theNDLDataA,Z which has A and Z of NDL data by T. Koi
32 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
33 // 101111 Add Special treatment for Be9(n,2n)Be8(2a) case by T. Koi
34 //
35 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
36 //
37 // June-2019 - E. Mendoza --> Added protection against residual with Z<0 or A<Z + adjust_final_state
38 // is not applied when data is in MF=6 format (no correlated particle emission) + bug correction
39 // (add Q value info to G4ParticleHPNBodyPhaseSpace).
40
41 #include "G4ParticleHPInelasticBaseFS.hh"
42
43 #include "G4Alpha.hh"
44 #include "G4Electron.hh"
45 #include "G4He3.hh"
46 #include "G4IonTable.hh"
47 #include "G4NucleiProperties.hh"
48 #include "G4Nucleus.hh"
49 #include "G4ParticleHPDataUsed.hh"
50 #include "G4ParticleHPManager.hh"
51
52 G4ParticleHPInelasticBaseFS::G4ParticleHPInelasticBaseFS()
53 {
54 hasXsec = true;
55 theXsection = new G4ParticleHPVector;
56 }
57
58 G4ParticleHPInelasticBaseFS::~G4ParticleHPInelasticBaseFS()
59 {
60 delete theXsection;
61 delete theEnergyDistribution;
62 delete theFinalStatePhotons;
63 delete theEnergyAngData;
64 delete theAngularDistribution;
65 }
66
67 void G4ParticleHPInelasticBaseFS::InitGammas(G4double AR, G4double ZR)
68 {
69 G4int Z = G4lrint(ZR);
70 G4int A = G4lrint(AR);
71 std::ostringstream ost;
72 ost << gammaPath << "z" << Z << ".a" << A;
73 G4String aName = ost.str();
74 std::ifstream from(aName, std::ios::in);
75
76 if (!from) return; // no data found for this isotope
77 std::ifstream theGammaData(aName, std::ios::in);
78
79 theNuclearMassDifference = G4NucleiProperties::GetBindingEnergy(A, Z)
80 - G4NucleiProperties::GetBindingEnergy(theBaseA, theBaseZ);
81 theGammas.Init(theGammaData);
82 }
83
84 void G4ParticleHPInelasticBaseFS::Init(G4double A, G4double Z, G4int M,
85 const G4String& dirName,
86 const G4String& bit, G4ParticleDefinition*)
87 {
88 gammaPath = fManager->GetNeutronHPPath() + "/Inelastic/Gammas/";
89 G4String tString = dirName;
90 SetA_Z(A, Z, M);
91 G4bool dbool = true;
92 G4ParticleHPDataUsed aFile =
93 theNames.GetName(theBaseA, theBaseZ, M, tString, bit, dbool);
94 SetAZMs(aFile);
95 G4String filename = aFile.GetName();
96 #ifdef G4VERBOSE
97 if (fManager->GetDEBUG())
98 G4cout << " G4ParticleHPInelasticBaseFS::Init FILE " << filename << G4endl;
99 #endif
100
101 if (!dbool || (theBaseZ <= 2 && (theNDLDataZ != theBaseZ || theNDLDataA != theBaseA)))
102 {
103 #ifdef G4PHPDEBUG
104 if (fManager->GetDEBUG())
105 G4cout << "Skipped = " << filename << " " << A << " " << Z << G4endl;
106 #endif
107 hasAnyData = false;
108 hasFSData = false;
109 hasXsec = false;
110 return;
111 }
112
113 std::istringstream theData(std::ios::in);
114 fManager->GetDataStream(filename, theData);
115
116 if (!theData) //"!" is a operator of ios
117 {
118 hasAnyData = false;
119 hasFSData = false;
120 hasXsec = false;
121 return; // no data for exactly this isotope and FS
122 }
123 // here we go
124 G4int infoType, dataType, dummy = INT_MAX;
125 hasFSData = false;
126 while (theData >> infoType) // Loop checking, 11.05.2015, T. Koi
127 {
128 theData >> dataType;
129
130 if (dummy == INT_MAX) theData >> Qvalue >> dummy;
131 Qvalue *= CLHEP::eV;
132 // In G4NDL4.5 this value is the MT number (<1000),
133 // in others is que Q-value in eV
134
135 if (dataType == 3) {
136 G4int total;
137 theData >> total;
138 theXsection->Init(theData, total, CLHEP::eV);
139 }
140 else if (dataType == 4) {
141 theAngularDistribution = new G4ParticleHPAngular;
142 theAngularDistribution->Init(theData);
143 hasFSData = true;
144 }
145 else if (dataType == 5) {
146 theEnergyDistribution = new G4ParticleHPEnergyDistribution;
147 theEnergyDistribution->Init(theData);
148 hasFSData = true;
149 }
150 else if (dataType == 6) {
151 theEnergyAngData = new G4ParticleHPEnAngCorrelation(theProjectile);
152 theEnergyAngData->Init(theData);
153 hasFSData = true;
154 }
155 else if (dataType == 12) {
156 theFinalStatePhotons = new G4ParticleHPPhotonDist;
157 theFinalStatePhotons->InitMean(theData);
158 hasFSData = true;
159 }
160 else if (dataType == 13) {
161 theFinalStatePhotons = new G4ParticleHPPhotonDist;
162 theFinalStatePhotons->InitPartials(theData, theXsection);
163 hasFSData = true;
164 }
165 else if (dataType == 14) {
166 theFinalStatePhotons->InitAngular(theData);
167 hasFSData = true;
168 }
169 else if (dataType == 15) {
170 theFinalStatePhotons->InitEnergies(theData);
171 hasFSData = true;
172 }
173 else {
174 G4ExceptionDescription ed;
175 ed << "Z=" << theBaseZ << " A=" << theBaseA << " dataType=" << dataType
176 << " projectile: " << theProjectile->GetParticleName();
177 G4Exception("G4ParticleHPInelasticBaseFS::Init", "hadr01", JustWarning,
178 ed, "Data-type unknown");
179 }
180 }
181 }
182
183 void G4ParticleHPInelasticBaseFS::BaseApply(const G4HadProjectile& theTrack,
184 G4ParticleDefinition** theDefs,
185 G4int nDef)
186 {
187 // G4cout << "G4ParticleHPInelasticBaseFS::BaseApply started" << G4endl;
188 // prepare neutron
189 if (theResult.Get() == nullptr) theResult.Put(new G4HadFinalState);
190 theResult.Get()->Clear();
191 G4double eKinetic = theTrack.GetKineticEnergy();
192 G4ReactionProduct incidReactionProduct(theTrack.GetDefinition());
193 incidReactionProduct.SetMomentum(theTrack.Get4Momentum().vect());
194 incidReactionProduct.SetKineticEnergy(eKinetic);
195
196 // prepare target
197 G4double targetMass = G4NucleiProperties::GetNuclearMass(theBaseA, theBaseZ)
198 /CLHEP::neutron_mass_c2;
199 /*
200 G4cout << " Target Z=" << theBaseZ << " A=" << theBaseA
201 << " projectile: " << theTrack.GetDefinition()->GetParticleName()
202 << " Ekin(MeV)=" << theTrack.GetKineticEnergy()
203 << " Mtar/Mn=" << targetMass
204 << " hasFS=" << HasFSData() << G4endl;
205 */
206 // give priority to ENDF vales for target mass (VI: to be understood!?)
207 if (theEnergyAngData != nullptr) {
208 G4double x = theEnergyAngData->GetTargetMass();
209 if(0.0 < x) { targetMass = x/CLHEP::neutron_mass_c2; }
210 }
211 if (theAngularDistribution != nullptr) {
212 G4double x = theAngularDistribution->GetTargetMass();
213 if(0.0 < x) { targetMass = x/CLHEP::neutron_mass_c2; }
214 }
215
216 // 110512 TKDB ENDF-VII.0 21Sc45 has trouble in MF4MT22 (n,np) targetMass is not properly
217 // recorded. TKDB targetMass = 0; ENDF-VII.0 21Sc45 has trouble in MF4MT22 (n,np)
218
219 G4Nucleus aNucleus;
220 G4ReactionProduct theTarget;
221
222 G4ThreeVector neuVelo =
223 incidReactionProduct.GetMomentum() / CLHEP::neutron_mass_c2;
224
225 theTarget =
226 aNucleus.GetBiasedThermalNucleus(targetMass, neuVelo,
227 theTrack.GetMaterial()->GetTemperature());
228 theTarget.SetDefinition(ionTable->GetIon(theBaseZ, theBaseA, 0.0));
229
230 // prepare energy in target rest frame
231 G4ReactionProduct boosted;
232 boosted.Lorentz(incidReactionProduct, theTarget);
233 eKinetic = boosted.GetKineticEnergy();
234
235 G4double orgMomentum = boosted.GetMomentum().mag();
236
237 // Take N-body phase-space distribution, if no other data present.
238 if (!HasFSData()) // adding the residual is trivial here @@@
239 {
240 G4ParticleHPNBodyPhaseSpace thePhaseSpaceDistribution;
241 G4double aPhaseMass = 0;
242 G4int ii;
243 for (ii = 0; ii < nDef; ++ii) {
244 aPhaseMass += theDefs[ii]->GetPDGMass();
245 }
246
247 //-----------------------------------------------------------------
248 if (std::abs(Qvalue) < CLHEP::keV) {
249 // Not in the G4NDL lib or not calculated yet:
250
251 // Calculate residual:
252 G4int ResidualA = theBaseA;
253 G4int ResidualZ = theBaseZ;
254 for (ii = 0; ii < nDef; ++ii) {
255 ResidualZ -= theDefs[ii]->GetAtomicNumber();
256 ResidualA -= theDefs[ii]->GetBaryonNumber();
257 }
258
259 if (ResidualA > 0 && ResidualZ > 0) {
260 G4ParticleDefinition* resid = ionTable->GetIon(ResidualZ, ResidualA);
261 Qvalue =
262 incidReactionProduct.GetMass() + theTarget.GetMass() - aPhaseMass - resid->GetPDGMass();
263 }
264
265 if (std::abs(Qvalue) > 400 * CLHEP::MeV) {
266 // Then Q value is probably too large ...
267 Qvalue = 1.1 * CLHEP::keV;
268 }
269 }
270 //----------------------------------------------------------------------------
271
272 thePhaseSpaceDistribution.Init(aPhaseMass, nDef);
273 thePhaseSpaceDistribution.SetProjectileRP(&incidReactionProduct);
274 thePhaseSpaceDistribution.SetTarget(&theTarget);
275 thePhaseSpaceDistribution.SetQValue(Qvalue);
276
277 for (ii = 0; ii < nDef; ++ii) {
278 G4double massCode = 1000. * std::abs(theDefs[ii]->GetPDGCharge());
279 massCode += theDefs[ii]->GetBaryonNumber();
280 G4double dummy = 0;
281 G4ReactionProduct* aSec = thePhaseSpaceDistribution.Sample(eKinetic, massCode, dummy);
282 aSec->Lorentz(*aSec, -1. * theTarget);
283 auto aPart = new G4DynamicParticle();
284 aPart->SetDefinition(aSec->GetDefinition());
285 aPart->SetMomentum(aSec->GetMomentum());
286 delete aSec;
287 theResult.Get()->AddSecondary(aPart, secID);
288 #ifdef G4VERBOSE
289 if (fManager->GetDEBUG())
290 G4cout << this << " G4ParticleHPInelasticBaseFS::BaseApply NoFSData add secondary "
291 << aPart->GetParticleDefinition()->GetParticleName()
292 << " E= " << aPart->GetKineticEnergy() << " NSECO "
293 << theResult.Get()->GetNumberOfSecondaries() << G4endl;
294 #endif
295 }
296
297 theResult.Get()->SetStatusChange(stopAndKill);
298 // Final momentum check should be done before return
299 const G4ParticleDefinition* targ_pd = ionTable->GetIon(theBaseZ, theBaseA, 0.0);
300 G4double mass = targ_pd->GetPDGMass();
301 G4LorentzVector targ_4p_lab(theTarget.GetMomentum(),
302 std::sqrt(mass*mass + theTarget.GetMomentum().mag2()));
303 G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
304 G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
305 adjust_final_state(init_4p_lab);
306 return;
307 }
308
309 // set target and neutron in the relevant exit channel
310 if (theAngularDistribution != nullptr) {
311 theAngularDistribution->SetTarget(theTarget);
312 theAngularDistribution->SetProjectileRP(incidReactionProduct);
313 }
314 else if (theEnergyAngData != nullptr) {
315 theEnergyAngData->SetTarget(theTarget);
316 theEnergyAngData->SetProjectileRP(incidReactionProduct);
317 }
318
319 G4ReactionProductVector* tmpHadrons = nullptr;
320 G4int dummy;
321 std::size_t i;
322
323 if (theEnergyAngData != nullptr) {
324 tmpHadrons = theEnergyAngData->Sample(eKinetic);
325 auto hproj = theTrack.GetDefinition();
326 G4int projA = hproj->GetBaryonNumber();
327 G4int projZ = G4lrint(hproj->GetPDGCharge()/CLHEP::eplus);
328
329 if (!fManager->GetDoNotAdjustFinalState()) {
330 // Adjust A and Z in the case of miss much between selected data and target nucleus
331 if (tmpHadrons != nullptr) {
332 G4int sumA = 0;
333 G4int sumZ = 0;
334 G4int maxA = 0;
335 G4int jAtMaxA = 0;
336 for (G4int j = 0; j != (G4int)tmpHadrons->size(); ++j) {
337 auto had = ((*tmpHadrons)[j])->GetDefinition();
338 if (had->GetBaryonNumber() > maxA) {
339 maxA = had->GetBaryonNumber();
340 jAtMaxA = j;
341 }
342 sumA += had->GetBaryonNumber();
343 sumZ += G4lrint(had->GetPDGCharge()/CLHEP::eplus);
344 }
345 G4int dA = theBaseA + projA - sumA;
346 G4int dZ = theBaseZ + projZ - sumZ;
347 if (dA < 0 || dZ < 0) {
348 auto p = ((*tmpHadrons)[jAtMaxA])->GetDefinition();
349 G4int newA = p->GetBaryonNumber() + dA;
350 G4int newZ = G4lrint(p->GetPDGCharge()/CLHEP::eplus) + dZ;
351 if (newA > newZ && newZ > 0) {
352 G4ParticleDefinition* pd = ionTable->GetIon(newZ, newA);
353 ((*tmpHadrons)[jAtMaxA])->SetDefinition(pd);
354 }
355 }
356 }
357 }
358 }
359 else if (theAngularDistribution != nullptr) {
360
361 auto Done = new G4bool[nDef];
362 G4int i0;
363 for (i0 = 0; i0 < nDef; ++i0)
364 Done[i0] = false;
365
366 tmpHadrons = new G4ReactionProductVector;
367 G4ReactionProduct* aHadron;
368 G4double localMass = G4NucleiProperties::GetNuclearMass(theBaseA, theBaseZ);
369 G4ThreeVector bufferedDirection(0, 0, 0);
370 for (i0 = 0; i0 < nDef; ++i0) {
371 if (!Done[i0]) {
372 aHadron = new G4ReactionProduct;
373 if (theEnergyDistribution != nullptr) {
374 aHadron->SetDefinition(theDefs[i0]);
375 aHadron->SetKineticEnergy(theEnergyDistribution->Sample(eKinetic, dummy));
376 }
377 else if (nDef == 1) {
378 aHadron->SetDefinition(theDefs[i0]);
379 aHadron->SetKineticEnergy(eKinetic);
380 }
381 else if (nDef == 2) {
382 aHadron->SetDefinition(theDefs[i0]);
383 aHadron->SetKineticEnergy(50 * CLHEP::MeV);
384 }
385 else {
386 throw G4HadronicException(
387 __FILE__, __LINE__,
388 "No energy distribution to sample from in InelasticBaseFS::BaseApply");
389 }
390 theAngularDistribution->SampleAndUpdate(*aHadron);
391 if (theEnergyDistribution == nullptr && nDef == 2) {
392 if (i0 == 0) {
393 G4double mass1 = theDefs[0]->GetPDGMass();
394 G4double mass2 = theDefs[1]->GetPDGMass();
395 G4double massn = theProjectile->GetPDGMass();
396 G4int z1 = theBaseZ -
397 G4lrint((theDefs[0]->GetPDGCharge() + theDefs[1]->GetPDGCharge())/CLHEP::eplus);
398 G4int a1 = theBaseA - theDefs[0]->GetBaryonNumber() - theDefs[1]->GetBaryonNumber();
399 G4double concreteMass = G4NucleiProperties::GetNuclearMass(a1, z1);
400 G4double availableEnergy = eKinetic + massn + localMass - mass1 - mass2
401 - concreteMass;
402 // available kinetic energy in CMS (non relativistic)
403 G4double emin = availableEnergy + mass1 + mass2
404 - std::sqrt((mass1 + mass2) * (mass1 + mass2) + orgMomentum * orgMomentum);
405 G4double p1 = std::sqrt(2. * mass2 * emin);
406 bufferedDirection = p1 * aHadron->GetMomentum().unit();
407 #ifdef G4PHPDEBUG
408 if (fManager->GetDEBUG()) // @@@@@ verify the nucleon counting...
409 {
410 G4cout << "G4ParticleHPInelasticBaseFS " << z1 << " " << theBaseZ
411 << " " << a1 << " " << theBaseA << " " << availableEnergy
412 << " " << emin << G4endl;
413 }
414 #endif
415 }
416 else {
417 bufferedDirection = -bufferedDirection;
418 }
419 // boost from cms to lab
420 aHadron->SetTotalEnergy(
421 std::sqrt(aHadron->GetMass() * aHadron->GetMass() + bufferedDirection.mag2()));
422 aHadron->SetMomentum(bufferedDirection);
423 aHadron->Lorentz(*aHadron, -1. * (theTarget + incidReactionProduct));
424 #ifdef G4PHPDEBUG
425 if (fManager->GetDEBUG()) {
426 G4cout << " G4ParticleHPInelasticBaseFS E=" << aHadron->GetTotalEnergy()
427 << " P=" << aHadron->GetMomentum()
428 << " bufDir^2=" << bufferedDirection.mag2()
429 << G4endl;
430 }
431 #endif
432 }
433 tmpHadrons->push_back(aHadron);
434 #ifdef G4PHPDEBUG
435 if (fManager->GetDEBUG())
436 G4cout << " G4ParticleHPInelasticBaseFS::BaseApply FSData add secondary "
437 << aHadron->GetDefinition()->GetParticleName()
438 << " E= " << aHadron->GetKineticEnergy() << G4endl;
439 #endif
440 }
441 }
442 delete[] Done;
443 }
444 else {
445 throw G4HadronicException(__FILE__, __LINE__, "No data to create InelasticFS");
446 }
447
448 G4ReactionProductVector* thePhotons = nullptr;
449 if (theFinalStatePhotons != nullptr) {
450 // the photon distributions are in the Nucleus rest frame.
451 G4ReactionProduct boosted_tmp;
452 boosted_tmp.Lorentz(incidReactionProduct, theTarget);
453 G4double anEnergy = boosted_tmp.GetKineticEnergy();
454 thePhotons = theFinalStatePhotons->GetPhotons(anEnergy);
455 if (thePhotons != nullptr) {
456 for (i = 0; i < thePhotons->size(); ++i) {
457 // back to lab
458 thePhotons->operator[](i)->Lorentz(*(thePhotons->operator[](i)), -1. * theTarget);
459 }
460 }
461 }
462 else if (theEnergyAngData != nullptr) {
463 // PA130927: do not create photons to adjust binding energy
464 G4bool bAdjustPhotons{true};
465 #ifdef PHP_AS_HP
466 bAdjustPhotons = true;
467 #else
468 if (fManager->GetDoNotAdjustFinalState()) bAdjustPhotons = false;
469 #endif
470
471 if (bAdjustPhotons) {
472 G4double theGammaEnergy = theEnergyAngData->GetTotalMeanEnergy();
473 G4double anEnergy = boosted.GetKineticEnergy();
474 theGammaEnergy = anEnergy - theGammaEnergy;
475 theGammaEnergy += theNuclearMassDifference;
476 G4double eBindProducts = 0;
477 G4double eBindN = 0;
478 G4double eBindP = 0;
479 G4double eBindD = G4NucleiProperties::GetBindingEnergy(2, 1);
480 G4double eBindT = G4NucleiProperties::GetBindingEnergy(3, 1);
481 G4double eBindHe3 = G4NucleiProperties::GetBindingEnergy(3, 2);
482 G4double eBindA = G4NucleiProperties::GetBindingEnergy(4, 2);
483 G4int ia = 0;
484 for (auto const & had : *tmpHadrons) {
485 if (had->GetDefinition() == G4Neutron::Neutron()) {
486 eBindProducts += eBindN;
487 }
488 else if (had->GetDefinition() == G4Proton::Proton()) {
489 eBindProducts += eBindP;
490 }
491 else if (had->GetDefinition() == G4Deuteron::Deuteron()) {
492 eBindProducts += eBindD;
493 }
494 else if (had->GetDefinition() == G4Triton::Triton()) {
495 eBindProducts += eBindT;
496 }
497 else if (had->GetDefinition() == G4He3::He3()) {
498 eBindProducts += eBindHe3;
499 }
500 else if (had->GetDefinition() == G4Alpha::Alpha()) {
501 eBindProducts += eBindA;
502 ia++;
503 }
504 }
505 theGammaEnergy += eBindProducts;
506
507 #ifdef G4PHPDEBUG
508 if (fManager->GetDEBUG())
509 G4cout << " G4ParticleHPInelasticBaseFS::BaseApply gamma Energy " << theGammaEnergy
510 << " eBindProducts " << eBindProducts << G4endl;
511 #endif
512
513 // Special treatment for Be9 + n -> 2n + Be8 -> 2n + a + a
514 if (theBaseZ == 4 && theBaseA == 9) {
515 // This only valid for G4NDL3.13,,,
516 if (ia == 2 && std::abs(G4NucleiProperties::GetBindingEnergy(8, 4) -
517 G4NucleiProperties::GetBindingEnergy(9, 4) -
518 theNuclearMassDifference) < CLHEP::keV)
519 {
520 theGammaEnergy -= (2 * eBindA);
521 }
522 }
523
524 if (theGammaEnergy > 0.0) {
525 for (G4int iLevel = theGammas.GetNumberOfLevels() - 1; iLevel > 0; --iLevel) {
526 G4double e = theGammas.GetLevelEnergy(iLevel);
527 if (e < theGammaEnergy) {
528 thePhotons = theGammas.GetDecayGammas(iLevel);
529 theGammaEnergy -= e;
530 break;
531 }
532 }
533 }
534 }
535 }
536 // fill the result
537 std::size_t nSecondaries = tmpHadrons->size();
538 std::size_t nPhotons = 0;
539 if (thePhotons != nullptr) {
540 nPhotons = thePhotons->size();
541 }
542 nSecondaries += nPhotons;
543
544 G4DynamicParticle* theSec;
545 #ifdef G4PHPDEBUG
546 if (fManager->GetDEBUG())
547 G4cout << " G4ParticleHPInelasticBaseFS::BaseApply N hadrons " << nSecondaries - nPhotons
548 << G4endl;
549 #endif
550
551 for (i = 0; i < nSecondaries - nPhotons; ++i) {
552 theSec = new G4DynamicParticle;
553 theSec->SetDefinition(tmpHadrons->operator[](i)->GetDefinition());
554 theSec->SetMomentum(tmpHadrons->operator[](i)->GetMomentum());
555 theResult.Get()->AddSecondary(theSec, secID);
556 #ifdef G4PHPDEBUG
557 if (fManager->GetDEBUG())
558 G4cout << " G4ParticleHPInelasticBaseFS::BaseApply add secondary2 "
559 << theSec->GetParticleDefinition()->GetParticleName()
560 << " E=" << theSec->GetKineticEnergy() << " Nsec="
561 << theResult.Get()->GetNumberOfSecondaries() << G4endl;
562 #endif
563 delete (*tmpHadrons)[i];
564 }
565 #ifdef G4PHPDEBUG
566 if (fManager->GetDEBUG())
567 G4cout << " G4ParticleHPInelasticBaseFS::BaseApply N photons " << nPhotons << G4endl;
568 #endif
569 if (thePhotons != nullptr) {
570 for (i = 0; i < nPhotons; ++i) {
571 theSec = new G4DynamicParticle;
572 theSec->SetDefinition((*thePhotons)[i]->GetDefinition());
573 theSec->SetMomentum((*thePhotons)[i]->GetMomentum());
574 theResult.Get()->AddSecondary(theSec, secID);
575 #ifdef G4PHPDEBUG
576 if (fManager->GetDEBUG())
577 G4cout << " G4ParticleHPInelasticBaseFS::BaseApply add secondary3 "
578 << theSec->GetParticleDefinition()->GetParticleName()
579 << " E=" << theSec->GetKineticEnergy() << " Nsec="
580 << theResult.Get()->GetNumberOfSecondaries() << G4endl;
581 #endif
582 delete thePhotons->operator[](i);
583 }
584 }
585
586 // some garbage collection
587 delete thePhotons;
588 delete tmpHadrons;
589
590 G4ParticleDefinition* targ_pd = ionTable->GetIon(theBaseZ, theBaseA, 0.0);
591 G4LorentzVector targ_4p_lab(
592 theTarget.GetMomentum(),
593 std::sqrt(targ_pd->GetPDGMass() * targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2()));
594 G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
595 G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
596
597 // if data in MF=6 format (no correlated particle emission), then adjust_final_state can give
598 // severe errors:
599 if (theEnergyAngData == nullptr) {
600 adjust_final_state(init_4p_lab);
601 }
602
603 // clean up the primary neutron
604 theResult.Get()->SetStatusChange(stopAndKill);
605 }
606