| 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 // 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi)
31 // 070606 bug fix and migrate to enable to Partial cases by T. Koi
32 // 080603 bug fix for Hadron Hyper News #932 by T. Koi
33 // 080612 bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4,6
34 // 080717 bug fix of calculation of residual momentum by T. Koi
35 // 080801 protect negative available energy by T. Koi
36 // introduce theNDLDataA,Z which has A and Z of NDL data by T. Koi
37 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
38 // 090514 Fix bug in IC electron emission case
39 // Contribution from Chao Zhang (Chao.Zhang@usd.edu) and Dongming Mei(Dongming.Mei@usd.edu)
40 // 100406 "nothingWasKnownOnHadron=1" then sample mu isotropic in CM
41 // add two_body_reaction
42 // 100909 add safty
43 // 101111 add safty for _nat_ data case in Binary reaction, but break conservation
44 // 110430 add Reaction Q value and break up flag (MF3::QI and LR)
45 //
46 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
47 // June-2019 - E. Mendoza - re-build "two_body_reaction", to be used by incident charged particles
48 // (now isotropic emission in the CMS). Also restrict nresp use below 20 MeV (for future
49 // developments). Add photon emission when no data available.
50 //
51 // nresp71_m03.hh and nresp71_m02.hh are alike. The only difference between m02 and m03
52 // is in the total carbon cross section that is properly included in the latter.
53 // These data are not used in nresp71_m0*.hh.
54 //
55
56 #include "G4ParticleHPInelasticCompFS.hh"
57
58 #include "G4Alpha.hh"
59 #include "G4Electron.hh"
60 #include "G4He3.hh"
61 #include "G4IonTable.hh"
62 #include "G4NRESP71M03.hh"
63 #include "G4NucleiProperties.hh"
64 #include "G4Nucleus.hh"
65 #include "G4ParticleHPDataUsed.hh"
66 #include "G4ParticleHPManager.hh"
67 #include "G4RandomDirection.hh"
68 #include "G4SystemOfUnits.hh"
69
70 #include <fstream>
71
72 G4ParticleHPInelasticCompFS::G4ParticleHPInelasticCompFS()
73 : G4ParticleHPFinalState()
74 {
75 QI.resize(51);
76 LR.resize(51);
77 for (G4int i = 0; i < 51; ++i) {
78 hasXsec = true;
79 theXsection[i] = nullptr;
80 theEnergyDistribution[i] = nullptr;
81 theAngularDistribution[i] = nullptr;
82 theEnergyAngData[i] = nullptr;
83 theFinalStatePhotons[i] = nullptr;
84 QI[i] = 0.0;
85 LR[i] = 0;
86 }
87 }
88
89 G4ParticleHPInelasticCompFS::~G4ParticleHPInelasticCompFS()
90 {
91 for (G4int i = 0; i < 51; ++i) {
92 if (theXsection[i] != nullptr) delete theXsection[i];
93 if (theEnergyDistribution[i] != nullptr) delete theEnergyDistribution[i];
94 if (theAngularDistribution[i] != nullptr) delete theAngularDistribution[i];
95 if (theEnergyAngData[i] != nullptr) delete theEnergyAngData[i];
96 if (theFinalStatePhotons[i] != nullptr) delete theFinalStatePhotons[i];
97 }
98 }
99
100 void G4ParticleHPInelasticCompFS::InitDistributionInitialState(G4ReactionProduct& inPart,
101 G4ReactionProduct& aTarget, G4int it)
102 {
103 if (theAngularDistribution[it] != nullptr) {
104 theAngularDistribution[it]->SetTarget(aTarget);
105 theAngularDistribution[it]->SetProjectileRP(inPart);
106 }
107
108 if (theEnergyAngData[it] != nullptr) {
109 theEnergyAngData[it]->SetTarget(aTarget);
110 theEnergyAngData[it]->SetProjectileRP(inPart);
111 }
112 }
113
114 void G4ParticleHPInelasticCompFS::InitGammas(G4double AR, G4double ZR)
115 {
116 G4int Z = G4lrint(ZR);
117 G4int A = G4lrint(AR);
118 std::ostringstream ost;
119 ost << gammaPath << "z" << Z << ".a" << A;
120 G4String aName = ost.str();
121 std::ifstream from(aName, std::ios::in);
122
123 if (!from) return; // no data found for this isotope
124 std::ifstream theGammaData(aName, std::ios::in);
125
126 theGammas.Init(theGammaData);
127 }
128
129 void G4ParticleHPInelasticCompFS::Init(G4double A, G4double Z, G4int M, const G4String& dirName,
130 const G4String& aFSType, G4ParticleDefinition*)
131 {
132 gammaPath = fManager->GetNeutronHPPath() + "/Inelastic/Gammas/";
133 const G4String& tString = dirName;
134 SetA_Z(A, Z, M);
135 G4bool dbool;
136 const G4ParticleHPDataUsed& aFile =
137 theNames.GetName(theBaseA, theBaseZ, M, tString, aFSType, dbool);
138 SetAZMs(aFile);
139 const G4String& filename = aFile.GetName();
140 #ifdef G4VERBOSE
141 if (fManager->GetDEBUG())
142 G4cout << " G4ParticleHPInelasticCompFS::Init FILE " << filename << G4endl;
143 #endif
144
145 SetAZMs(A, Z, M, aFile);
146
147 if (!dbool || (theBaseZ <= 2 && (theNDLDataZ != theBaseZ || theNDLDataA != theBaseA)))
148 {
149 #ifdef G4VERBOSE
150 if (fManager->GetDEBUG())
151 G4cout << "Skipped = " << filename << " " << A << " " << Z << G4endl;
152 #endif
153 hasAnyData = false;
154 hasFSData = false;
155 hasXsec = false;
156 return;
157 }
158 std::istringstream theData(std::ios::in);
159 fManager->GetDataStream(filename, theData);
160 if (!theData) //"!" is a operator of ios
161 {
162 hasAnyData = false;
163 hasFSData = false;
164 hasXsec = false;
165 return;
166 }
167 // here we go
168 G4int infoType, dataType, dummy;
169 G4int sfType, it;
170 hasFSData = false;
171 while (theData >> infoType) // Loop checking, 11.05.2015, T. Koi
172 {
173 hasFSData = true;
174 theData >> dataType;
175 theData >> sfType >> dummy;
176 it = 50;
177 if (sfType >= 600 || (sfType < 100 && sfType >= 50))
178 it = sfType % 50;
179 if (dataType == 3)
180 {
181 G4double dqi;
182 G4int ilr;
183 theData >> dqi >> ilr;
184
185 QI[it] = dqi * CLHEP::eV;
186 LR[it] = ilr;
187 theXsection[it] = new G4ParticleHPVector;
188 G4int total;
189 theData >> total;
190 theXsection[it]->Init(theData, total, CLHEP::eV);
191 }
192 else if (dataType == 4) {
193 theAngularDistribution[it] = new G4ParticleHPAngular;
194 theAngularDistribution[it]->Init(theData);
195 }
196 else if (dataType == 5) {
197 theEnergyDistribution[it] = new G4ParticleHPEnergyDistribution;
198 theEnergyDistribution[it]->Init(theData);
199 }
200 else if (dataType == 6) {
201 theEnergyAngData[it] = new G4ParticleHPEnAngCorrelation(theProjectile);
202 // G4cout << this << " CompFS theEnergyAngData " << it << theEnergyAngData[it] << G4endl;
203 theEnergyAngData[it]->Init(theData);
204 }
205 else if (dataType == 12) {
206 theFinalStatePhotons[it] = new G4ParticleHPPhotonDist;
207 theFinalStatePhotons[it]->InitMean(theData);
208 }
209 else if (dataType == 13) {
210 theFinalStatePhotons[it] = new G4ParticleHPPhotonDist;
211 theFinalStatePhotons[it]->InitPartials(theData, theXsection[50]);
212 }
213 else if (dataType == 14) {
214 theFinalStatePhotons[it]->InitAngular(theData);
215 }
216 else if (dataType == 15) {
217 theFinalStatePhotons[it]->InitEnergies(theData);
218 }
219 else {
220 G4ExceptionDescription ed;
221 ed << "Z=" << theBaseZ << " A=" << theBaseA << " dataType=" << dataType
222 << " projectile: " << theProjectile->GetParticleName();
223 G4Exception("G4ParticleHPInelasticCompFS::Init", "hadr01", JustWarning,
224 ed, "Data-type unknown");
225 }
226 }
227 }
228
229 G4int G4ParticleHPInelasticCompFS::SelectExitChannel(G4double eKinetic)
230 {
231 G4double running[50];
232 running[0] = 0;
233 G4int i;
234 for (i = 0; i < 50; ++i) {
235 if (i != 0) running[i] = running[i - 1];
236 if (theXsection[i] != nullptr) {
237 running[i] += std::max(0., theXsection[i]->GetXsec(eKinetic));
238 }
239 }
240 G4double random = G4UniformRand();
241 G4double sum = running[49];
242 G4int it = 50;
243 if (0 != sum) {
244 G4int i0;
245 for (i0 = 0; i0 < 50; ++i0) {
246 it = i0;
247 if (random < running[i0] / sum) break;
248 }
249 }
250 return it;
251 }
252
253 // n,p,d,t,he3,a
254 void G4ParticleHPInelasticCompFS::CompositeApply(const G4HadProjectile& theTrack,
255 G4ParticleDefinition* aDefinition)
256 {
257 // prepare neutron
258 if (theResult.Get() == nullptr) theResult.Put(new G4HadFinalState);
259 theResult.Get()->Clear();
260 G4double eKinetic = theTrack.GetKineticEnergy();
261 const G4HadProjectile* hadProjectile = &theTrack;
262 G4ReactionProduct incidReactionProduct(hadProjectile->GetDefinition());
263 incidReactionProduct.SetMomentum(hadProjectile->Get4Momentum().vect());
264 incidReactionProduct.SetKineticEnergy(eKinetic);
265
266 // prepare target
267 for (G4int i = 0; i < 50; ++i) {
268 if (theXsection[i] != nullptr) {
269 break;
270 }
271 }
272
273 G4double targetMass = G4NucleiProperties::GetNuclearMass(theBaseA, theBaseZ);
274 #ifdef G4VERBOSE
275 if (fManager->GetDEBUG())
276 G4cout << "G4ParticleHPInelasticCompFS::CompositeApply A=" << theBaseA << " Z="
277 << theBaseZ << " incident " << hadProjectile->GetDefinition()->GetParticleName()
278 << G4endl;
279 #endif
280 G4ReactionProduct theTarget;
281 G4Nucleus aNucleus;
282 // G4ThreeVector neuVelo =
283 // (1./hadProjectile->GetDefinition()->GetPDGMass())*incidReactionProduct.GetMomentum(); theTarget
284 // = aNucleus.GetBiasedThermalNucleus( targetMass/hadProjectile->GetDefinition()->GetPDGMass() ,
285 // neuVelo, theTrack.GetMaterial()->GetTemperature()); G4Nucleus::GetBiasedThermalNucleus requests
286 // normalization of mass and velocity in neutron mass
287 G4ThreeVector neuVelo = incidReactionProduct.GetMomentum() / CLHEP::neutron_mass_c2;
288 theTarget = aNucleus.GetBiasedThermalNucleus(targetMass / CLHEP::neutron_mass_c2,
289 neuVelo, theTrack.GetMaterial()->GetTemperature());
290
291 theTarget.SetDefinition(G4IonTable::GetIonTable()->GetIon(theBaseZ, theBaseA, 0.0));
292
293 // prepare the residual mass
294 G4double residualMass = 0;
295 G4int residualZ = theBaseZ +
296 G4lrint((theProjectile->GetPDGCharge() - aDefinition->GetPDGCharge())/CLHEP::eplus);
297 G4int residualA = theBaseA + theProjectile->GetBaryonNumber() - aDefinition->GetBaryonNumber();
298 residualMass = G4NucleiProperties::GetNuclearMass(residualA, residualZ);
299
300 // prepare energy in target rest frame
301 G4ReactionProduct boosted;
302 boosted.Lorentz(incidReactionProduct, theTarget);
303 eKinetic = boosted.GetKineticEnergy();
304
305 // select exit channel for composite FS class.
306 G4int it = SelectExitChannel(eKinetic);
307
308 // E. Mendoza (2018) -- to use JENDL/AN-2005
309 if (theEnergyDistribution[it] == nullptr && theAngularDistribution[it] == nullptr
310 && theEnergyAngData[it] == nullptr)
311 {
312 if (theEnergyDistribution[50] != nullptr || theAngularDistribution[50] != nullptr
313 || theEnergyAngData[50] != nullptr)
314 {
315 it = 50;
316 }
317 }
318
319 // set target and neutron in the relevant exit channel
320 InitDistributionInitialState(incidReactionProduct, theTarget, it);
321
322 //---------------------------------------------------------------------//
323 // Hook for NRESP71MODEL
324 if (fManager->GetUseNRESP71Model() && eKinetic < 20 * CLHEP::MeV) {
325 if (theBaseZ == 6) // If the reaction is with Carbon...
326 {
327 if (theProjectile == G4Neutron::Definition()) {
328 if (use_nresp71_model(aDefinition, it, theTarget, boosted)) return;
329 }
330 }
331 }
332 //---------------------------------------------------------------------//
333
334 G4ReactionProductVector* thePhotons = nullptr;
335 G4ReactionProductVector* theParticles = nullptr;
336 G4ReactionProduct aHadron;
337 aHadron.SetDefinition(aDefinition); // what if only cross-sections exist ==> Na 23 11 @@@@
338 G4double availableEnergy = incidReactionProduct.GetKineticEnergy()
339 + incidReactionProduct.GetMass() - aHadron.GetMass()
340 + (targetMass - residualMass);
341
342 if (availableEnergy < 0) {
343 availableEnergy = 0;
344 }
345 G4int nothingWasKnownOnHadron = 0;
346 G4double eGamm = 0;
347 G4int iLevel = -1;
348 // max gamma energy and index
349 G4int imaxEx = theGammas.GetNumberOfLevels() - 1;
350
351 // without photon has it = 0
352 if (50 == it) {
353 // Excitation level is not determined
354 aHadron.SetKineticEnergy(availableEnergy * residualMass / (aHadron.GetMass() + residualMass));
355
356 // TK add safty 100909
357 G4double p2 =
358 (aHadron.GetTotalEnergy() * aHadron.GetTotalEnergy() - aHadron.GetMass() * aHadron.GetMass());
359 G4double p = (p2 > 0.0) ? std::sqrt(p2) : 0.0;
360 aHadron.SetMomentum(p * incidReactionProduct.GetMomentum() /
361 incidReactionProduct.GetTotalMomentum());
362 }
363 else {
364 iLevel = imaxEx;
365 }
366
367 if (theAngularDistribution[it] != nullptr) // MF4
368 {
369 if (theEnergyDistribution[it] != nullptr) // MF5
370 {
371 //************************************************************
372 /*
373 aHadron.SetKineticEnergy(theEnergyDistribution[it]->Sample(eKinetic, dummy));
374 G4double eSecN = aHadron.GetKineticEnergy();
375 */
376 //************************************************************
377 // EMendoza --> maximum allowable energy should be taken into account.
378 G4double dqi = 0.0;
379 if (QI[it] < 0 || 849 < QI[it])
380 dqi = QI[it]; // For backword compatibility QI introduced since G4NDL3.15
381 G4double MaxEne = eKinetic + dqi;
382 G4double eSecN = 0.;
383
384 G4int icounter = 0;
385 G4int icounter_max = 1024;
386 G4int dummy = 0;
387 do {
388 ++icounter;
389 if (icounter > icounter_max) {
390 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of "
391 << __FILE__ << "." << G4endl;
392 break;
393 }
394 eSecN = theEnergyDistribution[it]->Sample(eKinetic, dummy);
395 } while (eSecN > MaxEne); // Loop checking, 11.05.2015, T. Koi
396 aHadron.SetKineticEnergy(eSecN);
397 //************************************************************
398 eGamm = eKinetic - eSecN;
399 for (iLevel = imaxEx; iLevel >= 0; --iLevel) {
400 if (theGammas.GetLevelEnergy(iLevel) < eGamm) break;
401 }
402 if (iLevel < imaxEx && iLevel >= 0) {
403 if (G4UniformRand() > 0.5) {
404 ++iLevel;
405 }
406 }
407 }
408 else {
409 G4double eExcitation = 0;
410 for (iLevel = imaxEx; iLevel >= 0; --iLevel) {
411 if (theGammas.GetLevelEnergy(iLevel) < eKinetic) break;
412 }
413
414 // Use QI value for calculating excitation energy of residual.
415 G4bool useQI = false;
416 G4double dqi = QI[it];
417 if (dqi < 0 || 849 < dqi) useQI = true; // Former libraries do not have values in this range
418
419 if (useQI) {
420 eExcitation = std::max(0., QI[0] - QI[it]); // Bug fix 2333
421
422 // Re-evaluate iLevel based on this eExcitation
423 iLevel = 0;
424 G4bool find = false;
425 const G4double level_tolerance = 1.0 * CLHEP::keV;
426
427 // VI: the first level is ground
428 if (0 < imaxEx) {
429 for (iLevel = 1; iLevel <= imaxEx; ++iLevel) {
430 G4double elevel = theGammas.GetLevelEnergy(iLevel);
431 if (std::abs(eExcitation - elevel) < level_tolerance) {
432 find = true;
433 break;
434 }
435 if (eExcitation < elevel) {
436 find = true;
437 iLevel = std::max(iLevel - 1, 0);
438 break;
439 }
440 }
441
442 // If proper level cannot be found, use the maximum level
443 if (!find) iLevel = imaxEx;
444 }
445 }
446
447 if (fManager->GetDEBUG() && eKinetic - eExcitation < 0) {
448 throw G4HadronicException(
449 __FILE__, __LINE__,
450 "SEVERE: InelasticCompFS: Consistency of data not good enough, please file report");
451 }
452 if (eKinetic - eExcitation < 0) eExcitation = 0;
453 if (iLevel != -1) aHadron.SetKineticEnergy(eKinetic - eExcitation);
454 }
455 theAngularDistribution[it]->SampleAndUpdate(aHadron);
456
457 if (theFinalStatePhotons[it] == nullptr) {
458 thePhotons = theGammas.GetDecayGammas(iLevel);
459 eGamm -= theGammas.GetLevelEnergy(iLevel);
460 }
461 }
462 else if (theEnergyAngData[it] != nullptr) // MF6
463 {
464 theParticles = theEnergyAngData[it]->Sample(eKinetic);
465
466 // Adjust A and Z in the case of miss much between selected data and target nucleus
467 if (theParticles != nullptr) {
468 G4int sumA = 0;
469 G4int sumZ = 0;
470 G4int maxA = 0;
471 G4int jAtMaxA = 0;
472 for (G4int j = 0; j != (G4int)theParticles->size(); ++j) {
473 auto ptr = theParticles->at(j);
474 G4int barnum = ptr->GetDefinition()->GetBaryonNumber();
475 if (barnum > maxA) {
476 maxA = barnum;
477 jAtMaxA = j;
478 }
479 sumA += barnum;
480 sumZ += G4lrint(ptr->GetDefinition()->GetPDGCharge()/CLHEP::eplus);
481 }
482 G4int dA = theBaseA + hadProjectile->GetDefinition()->GetBaryonNumber() - sumA;
483 G4int dZ = theBaseZ +
484 G4lrint(hadProjectile->GetDefinition()->GetPDGCharge()/CLHEP::eplus) - sumZ;
485 if (dA < 0 || dZ < 0) {
486 G4int newA = theParticles->at(jAtMaxA)->GetDefinition()->GetBaryonNumber() + dA;
487 G4int newZ =
488 G4lrint(theParticles->at(jAtMaxA)->GetDefinition()->GetPDGCharge()/CLHEP::eplus) + dZ;
489 G4ParticleDefinition* pd = ionTable->GetIon(newZ, newA);
490 theParticles->at(jAtMaxA)->SetDefinition(pd);
491 }
492 }
493 }
494 else {
495 // @@@ what to do, if we have photon data, but no info on the hadron itself
496 nothingWasKnownOnHadron = 1;
497 }
498
499 if (theFinalStatePhotons[it] != nullptr) {
500 // the photon distributions are in the Nucleus rest frame.
501 // TK residual rest frame
502 G4ReactionProduct boosted_tmp;
503 boosted_tmp.Lorentz(incidReactionProduct, theTarget);
504 G4double anEnergy = boosted_tmp.GetKineticEnergy();
505 thePhotons = theFinalStatePhotons[it]->GetPhotons(anEnergy);
506 G4double aBaseEnergy = theFinalStatePhotons[it]->GetLevelEnergy();
507 G4double testEnergy = 0;
508 if (thePhotons != nullptr && !thePhotons->empty()) {
509 aBaseEnergy -= (*thePhotons)[0]->GetTotalEnergy();
510 }
511 if (theFinalStatePhotons[it]->NeedsCascade()) {
512 while (aBaseEnergy > 0.01 * CLHEP::keV) // Loop checking, 11.05.2015, T. Koi
513 {
514 // cascade down the levels
515 G4bool foundMatchingLevel = false;
516 G4int closest = 2;
517 G4double deltaEold = -1;
518 for (G4int j = 1; j < it; ++j) {
519 if (theFinalStatePhotons[j] != nullptr) {
520 testEnergy = theFinalStatePhotons[j]->GetLevelEnergy();
521 }
522 else {
523 testEnergy = 0;
524 }
525 G4double deltaE = std::abs(testEnergy - aBaseEnergy);
526 if (deltaE < 0.1 * CLHEP::keV) {
527 G4ReactionProductVector* theNext = theFinalStatePhotons[j]->GetPhotons(anEnergy);
528 if (thePhotons != nullptr) thePhotons->push_back(theNext->operator[](0));
529 aBaseEnergy = testEnergy - theNext->operator[](0)->GetTotalEnergy();
530 delete theNext;
531 foundMatchingLevel = true;
532 break; // ===>
533 }
534 if (theFinalStatePhotons[j] != nullptr && (deltaE < deltaEold || deltaEold < 0.)) {
535 closest = j;
536 deltaEold = deltaE;
537 }
538 } // <=== the break goes here.
539 if (!foundMatchingLevel) {
540 G4ReactionProductVector* theNext = theFinalStatePhotons[closest]->GetPhotons(anEnergy);
541 if (thePhotons != nullptr) thePhotons->push_back(theNext->operator[](0));
542 aBaseEnergy = aBaseEnergy - theNext->operator[](0)->GetTotalEnergy();
543 delete theNext;
544 }
545 }
546 }
547 }
548
549 if (thePhotons != nullptr) {
550 for (auto const & p : *thePhotons) {
551 // back to lab
552 p->Lorentz(*p, -1. * theTarget);
553 }
554 }
555 if (nothingWasKnownOnHadron != 0) {
556 // In this case, hadron should be isotropic in CM
557 // Next 12 lines are Emilio's replacement
558 // G4double QM=(incidReactionProduct.GetMass()+targetMass)-(aHadron.GetMass()+residualMass);
559 // G4double eExcitation = QM-QI[it];
560 // G4double eExcitation = QI[0] - QI[it]; // Fix of bug #1838
561 // if(eExcitation<20*CLHEP::keV){eExcitation=0;}
562
563 G4double eExcitation = std::max(0., QI[0] - QI[it]); // Fix of bug #2333
564
565 two_body_reaction(&incidReactionProduct, &theTarget, &aHadron, eExcitation);
566 if (thePhotons == nullptr && eExcitation > 0) {
567 for (iLevel = imaxEx; iLevel >= 0; --iLevel) {
568 if (theGammas.GetLevelEnergy(iLevel) < eExcitation + 5 * keV) break; // 5 keV tolerance
569 }
570 thePhotons = theGammas.GetDecayGammas(iLevel);
571 }
572 }
573
574 // fill the result
575 // Beware - the recoil is not necessarily in the particles...
576 // Can be calculated from momentum conservation?
577 // The idea is that the particles ar emitted forst, and the gammas only once the
578 // recoil is on the residual; assumption is that gammas do not contribute to
579 // the recoil.
580 // This needs more design @@@
581
582 G4bool needsSeparateRecoil = false;
583 G4int totalBaryonNumber = 0;
584 G4int totalCharge = 0;
585 G4ThreeVector totalMomentum(0);
586 if (theParticles != nullptr) {
587 const G4ParticleDefinition* aDef;
588 for (std::size_t ii0 = 0; ii0 < theParticles->size(); ++ii0) {
589 aDef = (*theParticles)[ii0]->GetDefinition();
590 totalBaryonNumber += aDef->GetBaryonNumber();
591 totalCharge += G4lrint(aDef->GetPDGCharge()/CLHEP::eplus);
592 totalMomentum += (*theParticles)[ii0]->GetMomentum();
593 }
594 if (totalBaryonNumber
595 != theBaseA + hadProjectile->GetDefinition()->GetBaryonNumber())
596 {
597 needsSeparateRecoil = true;
598 residualA = theBaseA + hadProjectile->GetDefinition()->GetBaryonNumber()
599 - totalBaryonNumber;
600 residualZ = theBaseZ +
601 G4lrint((hadProjectile->GetDefinition()->GetPDGCharge() - totalCharge)/CLHEP::eplus);
602 }
603 }
604
605 std::size_t nPhotons = 0;
606 if (thePhotons != nullptr) {
607 nPhotons = thePhotons->size();
608 }
609
610 G4DynamicParticle* theSec;
611
612 if (theParticles == nullptr) {
613 theSec = new G4DynamicParticle;
614 theSec->SetDefinition(aHadron.GetDefinition());
615 theSec->SetMomentum(aHadron.GetMomentum());
616 theResult.Get()->AddSecondary(theSec, secID);
617 #ifdef G4VERBOSE
618 if (fManager->GetDEBUG())
619 G4cout << " G4ParticleHPInelasticCompFS::BaseApply add secondary1 "
620 << theSec->GetParticleDefinition()->GetParticleName()
621 << " E= " << theSec->GetKineticEnergy() << " NSECO "
622 << theResult.Get()->GetNumberOfSecondaries() << G4endl;
623 #endif
624
625 aHadron.Lorentz(aHadron, theTarget);
626 G4ReactionProduct theResidual;
627 theResidual.SetDefinition(ionTable->GetIon(residualZ, residualA, 0));
628 theResidual.SetKineticEnergy(aHadron.GetKineticEnergy() * aHadron.GetMass()
629 / theResidual.GetMass());
630
631 // 080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #6
632 // theResidual.SetMomentum(-1.*aHadron.GetMomentum());
633 G4ThreeVector incidentNeutronMomentum = incidReactionProduct.GetMomentum();
634 theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
635
636 theResidual.Lorentz(theResidual, -1. * theTarget);
637 G4ThreeVector totalPhotonMomentum(0, 0, 0);
638 if (thePhotons != nullptr) {
639 for (std::size_t i = 0; i < nPhotons; ++i) {
640 totalPhotonMomentum += (*thePhotons)[i]->GetMomentum();
641 }
642 }
643 theSec = new G4DynamicParticle;
644 theSec->SetDefinition(theResidual.GetDefinition());
645 theSec->SetMomentum(theResidual.GetMomentum() - totalPhotonMomentum);
646 theResult.Get()->AddSecondary(theSec, secID);
647 #ifdef G4VERBOSE
648 if (fManager->GetDEBUG())
649 G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary2 "
650 << theSec->GetParticleDefinition()->GetParticleName()
651 << " E= " << theSec->GetKineticEnergy() << " NSECO "
652 << theResult.Get()->GetNumberOfSecondaries() << G4endl;
653 #endif
654 }
655 else {
656 for (std::size_t i0 = 0; i0 < theParticles->size(); ++i0) {
657 theSec = new G4DynamicParticle;
658 theSec->SetDefinition((*theParticles)[i0]->GetDefinition());
659 theSec->SetMomentum((*theParticles)[i0]->GetMomentum());
660 theResult.Get()->AddSecondary(theSec, secID);
661 #ifdef G4VERBOSE
662 if (fManager->GetDEBUG())
663 G4cout << " G4ParticleHPInelasticCompFS::BaseApply add secondary3 "
664 << theSec->GetParticleDefinition()->GetParticleName()
665 << " E= " << theSec->GetKineticEnergy() << " NSECO "
666 << theResult.Get()->GetNumberOfSecondaries() << G4endl;
667 #endif
668 delete (*theParticles)[i0];
669 }
670 delete theParticles;
671 if (needsSeparateRecoil && residualZ != 0) {
672 G4ReactionProduct theResidual;
673 theResidual.SetDefinition(ionTable->GetIon(residualZ, residualA, 0));
674 G4double resiualKineticEnergy = theResidual.GetMass() * theResidual.GetMass();
675 resiualKineticEnergy += totalMomentum * totalMomentum;
676 resiualKineticEnergy = std::sqrt(resiualKineticEnergy) - theResidual.GetMass();
677 theResidual.SetKineticEnergy(resiualKineticEnergy);
678
679 // 080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4
680 // theResidual.SetMomentum(-1.*totalMomentum);
681 // G4ThreeVector incidentNeutronMomentum = incidReactionProduct.GetMomentum();
682 // theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
683 // 080717 TK Comment still do NOT include photon's mometum which produce by thePhotons
684 theResidual.SetMomentum(incidReactionProduct.GetMomentum() + theTarget.GetMomentum()
685 - totalMomentum);
686
687 theSec = new G4DynamicParticle;
688 theSec->SetDefinition(theResidual.GetDefinition());
689 theSec->SetMomentum(theResidual.GetMomentum());
690 theResult.Get()->AddSecondary(theSec, secID);
691 #ifdef G4VERBOSE
692 if (fManager->GetDEBUG())
693 G4cout << " G4ParticleHPInelasticCompFS::BaseApply add secondary4 "
694 << theSec->GetParticleDefinition()->GetParticleName()
695 << " E= " << theSec->GetKineticEnergy() << " NSECO "
696 << theResult.Get()->GetNumberOfSecondaries() << G4endl;
697 #endif
698 }
699 }
700 if (thePhotons != nullptr) {
701 for (std::size_t i = 0; i < nPhotons; ++i) {
702 theSec = new G4DynamicParticle;
703 // Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25,
704 // 2009 theSec->SetDefinition(G4Gamma::Gamma());
705 theSec->SetDefinition((*thePhotons)[i]->GetDefinition());
706 // But never cause real effect at least with G4NDL3.13 TK
707 theSec->SetMomentum((*thePhotons)[i]->GetMomentum());
708 theResult.Get()->AddSecondary(theSec, secID);
709 #ifdef G4VERBOSE
710 if (fManager->GetDEBUG())
711 G4cout << " G4ParticleHPInelasticCompFS::BaseApply add secondary5 "
712 << theSec->GetParticleDefinition()->GetParticleName()
713 << " E= " << theSec->GetKineticEnergy() << " NSECO "
714 << theResult.Get()->GetNumberOfSecondaries() << G4endl;
715 #endif
716
717 delete thePhotons->operator[](i);
718 }
719 // some garbage collection
720 delete thePhotons;
721 }
722
723 G4ParticleDefinition* targ_pd = ionTable->GetIon(theBaseZ, theBaseA, 0.0);
724 G4LorentzVector targ_4p_lab(
725 theTarget.GetMomentum(),
726 std::sqrt(targ_pd->GetPDGMass() * targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2()));
727 G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
728 G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
729 adjust_final_state(init_4p_lab);
730
731 // clean up the primary neutron
732 theResult.Get()->SetStatusChange(stopAndKill);
733 }
734
735 // Re-implemented by E. Mendoza (2019). Isotropic emission in the CMS:
736 // proj: projectile in target-rest-frame (input)
737 // targ: target in target-rest-frame (input)
738 // product: secondary particle in target-rest-frame (output)
739 // resExcitationEnergy: excitation energy of the residual nucleus
740 //
741 void G4ParticleHPInelasticCompFS::two_body_reaction(G4ReactionProduct* proj,
742 G4ReactionProduct* targ,
743 G4ReactionProduct* product,
744 G4double resExcitationEnergy)
745 {
746 // CMS system:
747 G4ReactionProduct theCMS = *proj + *targ;
748
749 // Residual definition:
750 G4int resZ = G4lrint((proj->GetDefinition()->GetPDGCharge() + targ->GetDefinition()->GetPDGCharge()
751 - product->GetDefinition()->GetPDGCharge())/CLHEP::eplus);
752 G4int resA = proj->GetDefinition()->GetBaryonNumber() + targ->GetDefinition()->GetBaryonNumber()
753 - product->GetDefinition()->GetBaryonNumber();
754 G4ReactionProduct theResidual;
755 theResidual.SetDefinition(ionTable->GetIon(resZ, resA, 0.0));
756
757 // CMS system:
758 G4ReactionProduct theCMSproj;
759 G4ReactionProduct theCMStarg;
760 theCMSproj.Lorentz(*proj, theCMS);
761 theCMStarg.Lorentz(*targ, theCMS);
762 // final Momentum in the CMS:
763 G4double totE = std::sqrt(theCMSproj.GetMass() * theCMSproj.GetMass()
764 + theCMSproj.GetTotalMomentum() * theCMSproj.GetTotalMomentum())
765 + std::sqrt(theCMStarg.GetMass() * theCMStarg.GetMass()
766 + theCMStarg.GetTotalMomentum() * theCMStarg.GetTotalMomentum());
767 G4double prodmass = product->GetMass();
768 G4double resmass = theResidual.GetMass() + resExcitationEnergy;
769 G4double fmomsquared = (totE * totE - (prodmass - resmass) * (prodmass - resmass)) *
770 (totE * totE - (prodmass + resmass) * (prodmass + resmass)) / (4.*totE*totE);
771 G4double fmom = (fmomsquared > 0) ? std::sqrt(fmomsquared) : 0.0;
772
773 // random (isotropic direction):
774 product->SetMomentum(fmom * G4RandomDirection());
775 product->SetTotalEnergy(std::sqrt(prodmass * prodmass + fmom * fmom)); // CMS
776 // Back to the LAB system:
777 product->Lorentz(*product, -1. * theCMS);
778 }
779
780 G4bool G4ParticleHPInelasticCompFS::use_nresp71_model(const G4ParticleDefinition* aDefinition,
781 const G4int itt,
782 const G4ReactionProduct& theTarget,
783 G4ReactionProduct& boosted)
784 {
785 if (aDefinition == G4Neutron::Definition()) // If the outgoing particle is a neutron...
786 {
787 // LR: flag LR in ENDF. It indicates whether there is breakup of the residual nucleus or not.
788 // it: exit channel (index of the carbon excited state)
789
790 // Added by A. R. Garcia (CIEMAT) to include the physics of C(N,N'3A) reactions from NRESP71.
791
792 if (LR[itt] > 0) // If there is breakup of the residual nucleus LR(flag LR in ENDF)>0 (i.e. Z=6
793 // MT=52-91 (it=MT-50)).
794 {
795 // Defining carbon as the target in the reference frame at rest.
796 G4ReactionProduct theCarbon(theTarget);
797
798 theCarbon.SetMomentum(G4ThreeVector());
799 theCarbon.SetKineticEnergy(0.);
800
801 // Creating four reaction products.
802 G4ReactionProduct theProds[4];
803
804 // Applying C(N,N'3A) reaction mechanisms in the target rest frame.
805 if (itt == 41) {
806 // QI=QM=-7.275 MeV for C-0(N,N')C-C(3A) in ENDF/B-VII.1.
807 // This is not the value of the QI of the first step according
808 // to the model. So we don't take it. Instead, we set the one
809 // we have calculated: QI=(mn+m12C)-(ma+m9Be+Ex9Be)=-8.130 MeV.
810 nresp71_model.ApplyMechanismI_NBeA2A(boosted, theCarbon, theProds, -8.130 /*QI[it]*/);
811 // N+C --> A[0]+9BE* | 9BE* --> N[1]+8BE | 8BE --> 2*A[2,3].
812 }
813 else {
814 nresp71_model.ApplyMechanismII_ACN2A(boosted, theCarbon, theProds, QI[itt]);
815 // N+C --> N'[0]+C* | C* --> A[1]+8BE | 8BE --> 2*A[2,3].
816 }
817
818 // Returning to the reference frame where the target was in motion.
819 for (auto& theProd : theProds) {
820 theProd.Lorentz(theProd, -1. * theTarget);
821 theResult.Get()->AddSecondary(
822 new G4DynamicParticle(theProd.GetDefinition(), theProd.GetMomentum()), secID);
823 }
824
825 // Killing the primary neutron.
826 theResult.Get()->SetStatusChange(stopAndKill);
827
828 return true;
829 }
830 }
831 else if (aDefinition == G4Alpha::Definition()) // If the outgoing particle is an alpha, ...
832 {
833 // Added by A. R. Garcia (CIEMAT) to include the physics of C(N,A)9BE reactions from NRESP71.
834
835 if (LR[itt] == 0) // If Z=6, an alpha particle is emitted and there is no breakup of the
836 // residual nucleus LR(flag LR in ENDF)==0.
837 {
838 // Defining carbon as the target in the reference frame at rest.
839 G4ReactionProduct theCarbon(theTarget);
840 theCarbon.SetMomentum(G4ThreeVector());
841 theCarbon.SetKineticEnergy(0.);
842
843 // Creating four reaction products.
844 G4ReactionProduct theProds[2];
845
846 // Applying C(N,A)9BE reaction mechanism.
847 nresp71_model.ApplyMechanismABE(boosted, theCarbon, theProds);
848 // N+C --> A[0]+9BE[1].
849
850 for (auto& theProd : theProds) {
851 // Returning to the system of reference where the target was in motion.
852 theProd.Lorentz(theProd, -1. * theTarget);
853 theResult.Get()->AddSecondary(
854 new G4DynamicParticle(theProd.GetDefinition(), theProd.GetMomentum()), secID);
855 }
856
857 // Killing the primary neutron.
858 theResult.Get()->SetStatusChange(stopAndKill);
859 return true;
860 }
861 G4Exception("G4ParticleHPInelasticCompFS::CompositeApply()", "G4ParticleInelasticCompFS.cc",
862 FatalException, "Alpha production with LR!=0.");
863 }
864 return false;
865 }
866