| 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 // particle_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans
29 //
30 // 070523 bug fix for G4FPE_DEBUG on by A. How
31 // 070606 bug fix and migrate to enable to Par
32 // 080603 bug fix for Hadron Hyper News #932 b
33 // 080612 bug fix contribution from Benoit Pir
34 // 080717 bug fix of calculation of residual m
35 // 080801 protect negative available energy by
36 // introduce theNDLDataA,Z which has A
37 // 081024 G4NucleiPropertiesTable:: to G4Nucle
38 // 090514 Fix bug in IC electron emission case
39 // Contribution from Chao Zhang (Chao.Z
40 // 100406 "nothingWasKnownOnHadron=1" then sam
41 // add two_body_reaction
42 // 100909 add safty
43 // 101111 add safty for _nat_ data case in Bin
44 // 110430 add Reaction Q value and break up fl
45 //
46 // P. Arce, June-2014 Conversion neutron_hp to
47 // June-2019 - E. Mendoza - re-build "two_body
48 // (now isotropic emission in the CMS). Also r
49 // developments). Add photon emission when no
50 //
51 // nresp71_m03.hh and nresp71_m02.hh are alike
52 // is in the total carbon cross section that i
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::G4ParticleHPInela
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::~G4ParticleHPInel
90 {
91 for (G4int i = 0; i < 51; ++i) {
92 if (theXsection[i] != nullptr) delete theX
93 if (theEnergyDistribution[i] != nullptr) d
94 if (theAngularDistribution[i] != nullptr)
95 if (theEnergyAngData[i] != nullptr) delete
96 if (theFinalStatePhotons[i] != nullptr) de
97 }
98 }
99
100 void G4ParticleHPInelasticCompFS::InitDistribu
101 G4ReactionPr
102 {
103 if (theAngularDistribution[it] != nullptr) {
104 theAngularDistribution[it]->SetTarget(aTar
105 theAngularDistribution[it]->SetProjectileR
106 }
107
108 if (theEnergyAngData[it] != nullptr) {
109 theEnergyAngData[it]->SetTarget(aTarget);
110 theEnergyAngData[it]->SetProjectileRP(inPa
111 }
112 }
113
114 void G4ParticleHPInelasticCompFS::InitGammas(G
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 thi
124 std::ifstream theGammaData(aName, std::ios::
125
126 theGammas.Init(theGammaData);
127 }
128
129 void G4ParticleHPInelasticCompFS::Init(G4doubl
130 const G
131 {
132 gammaPath = fManager->GetNeutronHPPath() + "
133 const G4String& tString = dirName;
134 SetA_Z(A, Z, M);
135 G4bool dbool;
136 const G4ParticleHPDataUsed& aFile =
137 theNames.GetName(theBaseA, theBaseZ, M, tS
138 SetAZMs(aFile);
139 const G4String& filename = aFile.GetName();
140 #ifdef G4VERBOSE
141 if (fManager->GetDEBUG())
142 G4cout << " G4ParticleHPInelasticCompFS::I
143 #endif
144
145 SetAZMs(A, Z, M, aFile);
146
147 if (!dbool || (theBaseZ <= 2 && (theNDLDataZ
148 {
149 #ifdef G4VERBOSE
150 if (fManager->GetDEBUG())
151 G4cout << "Skipped = " << filename << "
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 checkin
172 {
173 hasFSData = true;
174 theData >> dataType;
175 theData >> sfType >> dummy;
176 it = 50;
177 if (sfType >= 600 || (sfType < 100 && sfTy
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, CL
191 }
192 else if (dataType == 4) {
193 theAngularDistribution[it] = new G4Parti
194 theAngularDistribution[it]->Init(theData
195 }
196 else if (dataType == 5) {
197 theEnergyDistribution[it] = new G4Partic
198 theEnergyDistribution[it]->Init(theData)
199 }
200 else if (dataType == 6) {
201 theEnergyAngData[it] = new G4ParticleHPE
202 // G4cout << this << " CompFS theEn
203 theEnergyAngData[it]->Init(theData);
204 }
205 else if (dataType == 12) {
206 theFinalStatePhotons[it] = new G4Particl
207 theFinalStatePhotons[it]->InitMean(theDa
208 }
209 else if (dataType == 13) {
210 theFinalStatePhotons[it] = new G4Particl
211 theFinalStatePhotons[it]->InitPartials(t
212 }
213 else if (dataType == 14) {
214 theFinalStatePhotons[it]->InitAngular(th
215 }
216 else if (dataType == 15) {
217 theFinalStatePhotons[it]->InitEnergies(t
218 }
219 else {
220 G4ExceptionDescription ed;
221 ed << "Z=" << theBaseZ << " A=" << theBa
222 << " projectile: " << theProjectile->GetPar
223 G4Exception("G4ParticleHPInelasticCompFS
224 ed, "Data-type unknown");
225 }
226 }
227 }
228
229 G4int G4ParticleHPInelasticCompFS::SelectExitC
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
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::CompositeApp
255
256 {
257 // prepare neutron
258 if (theResult.Get() == nullptr) theResult.Pu
259 theResult.Get()->Clear();
260 G4double eKinetic = theTrack.GetKineticEnerg
261 const G4HadProjectile* hadProjectile = &theT
262 G4ReactionProduct incidReactionProduct(hadPr
263 incidReactionProduct.SetMomentum(hadProjecti
264 incidReactionProduct.SetKineticEnergy(eKinet
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::Ge
274 #ifdef G4VERBOSE
275 if (fManager->GetDEBUG())
276 G4cout << "G4ParticleHPInelasticCompFS::Co
277 << theBaseZ << " incident " << hadProject
278 << G4endl;
279 #endif
280 G4ReactionProduct theTarget;
281 G4Nucleus aNucleus;
282 // G4ThreeVector neuVelo =
283 // (1./hadProjectile->GetDefinition()->GetPD
284 // = aNucleus.GetBiasedThermalNucleus( targe
285 // neuVelo, theTrack.GetMaterial()->GetTempe
286 // normalization of mass and velocity in neu
287 G4ThreeVector neuVelo = incidReactionProduct
288 theTarget = aNucleus.GetBiasedThermalNucleus
289
290
291 theTarget.SetDefinition(G4IonTable::GetIonTa
292
293 // prepare the residual mass
294 G4double residualMass = 0;
295 G4int residualZ = theBaseZ +
296 G4lrint((theProjectile->GetPDGCharge() - a
297 G4int residualA = theBaseA + theProjectile->
298 residualMass = G4NucleiProperties::GetNuclea
299
300 // prepare energy in target rest frame
301 G4ReactionProduct boosted;
302 boosted.Lorentz(incidReactionProduct, theTar
303 eKinetic = boosted.GetKineticEnergy();
304
305 // select exit channel for composite FS clas
306 G4int it = SelectExitChannel(eKinetic);
307
308 // E. Mendoza (2018) -- to use JENDL/AN-2005
309 if (theEnergyDistribution[it] == nullptr &&
310 && theEnergyAngData[it] == nullptr)
311 {
312 if (theEnergyDistribution[50] != nullptr |
313 || theEnergyAngData[50] != nullptr)
314 {
315 it = 50;
316 }
317 }
318
319 // set target and neutron in the relevant ex
320 InitDistributionInitialState(incidReactionPr
321
322 //------------------------------------------
323 // Hook for NRESP71MODEL
324 if (fManager->GetUseNRESP71Model() && eKinet
325 if (theBaseZ == 6) // If the reaction is
326 {
327 if (theProjectile == G4Neutron::Definiti
328 if (use_nresp71_model(aDefinition, it,
329 }
330 }
331 }
332 //------------------------------------------
333
334 G4ReactionProductVector* thePhotons = nullpt
335 G4ReactionProductVector* theParticles = null
336 G4ReactionProduct aHadron;
337 aHadron.SetDefinition(aDefinition); // what
338 G4double availableEnergy = incidReactionProd
339 + incidReactionPr
340 + (targetMass - r
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()
350
351 // without photon has it = 0
352 if (50 == it) {
353 // Excitation level is not determined
354 aHadron.SetKineticEnergy(availableEnergy *
355
356 // TK add safty 100909
357 G4double p2 =
358 (aHadron.GetTotalEnergy() * aHadron.GetT
359 G4double p = (p2 > 0.0) ? std::sqrt(p2) :
360 aHadron.SetMomentum(p * incidReactionProdu
361 incidReactionProduct.GetTotalMomentum())
362 }
363 else {
364 iLevel = imaxEx;
365 }
366
367 if (theAngularDistribution[it] != nullptr)
368 {
369 if (theEnergyDistribution[it] != nullptr)
370 {
371 //**************************************
372 /*
373 aHadron.SetKineticEnergy(theEnergy
374 G4double eSecN = aHadron.GetKineti
375 */
376 //**************************************
377 // EMendoza --> maximum allowable energy
378 G4double dqi = 0.0;
379 if (QI[it] < 0 || 849 < QI[it])
380 dqi = QI[it]; // For backword compati
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
391 << __FILE__ << "." << G4endl;
392 break;
393 }
394 eSecN = theEnergyDistribution[it]->Sam
395 } while (eSecN > MaxEne); // Loop check
396 aHadron.SetKineticEnergy(eSecN);
397 //**************************************
398 eGamm = eKinetic - eSecN;
399 for (iLevel = imaxEx; iLevel >= 0; --iLe
400 if (theGammas.GetLevelEnergy(iLevel) <
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; --iLe
411 if (theGammas.GetLevelEnergy(iLevel) <
412 }
413
414 // Use QI value for calculating excitati
415 G4bool useQI = false;
416 G4double dqi = QI[it];
417 if (dqi < 0 || 849 < dqi) useQI = true;
418
419 if (useQI) {
420 eExcitation = std::max(0., QI[0] - QI[
421
422 // Re-evaluate iLevel based on this eE
423 iLevel = 0;
424 G4bool find = false;
425 const G4double level_tolerance = 1.0 *
426
427 // VI: the first level is ground
428 if (0 < imaxEx) {
429 for (iLevel = 1; iLevel <= imaxEx; +
430 G4double elevel = theGammas.GetLev
431 if (std::abs(eExcitation - elevel)
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,
443 if (!find) iLevel = imaxEx;
444 }
445 }
446
447 if (fManager->GetDEBUG() && eKinetic - e
448 throw G4HadronicException(
449 __FILE__, __LINE__,
450 "SEVERE: InelasticCompFS: Consistenc
451 }
452 if (eKinetic - eExcitation < 0) eExcitat
453 if (iLevel != -1) aHadron.SetKineticEner
454 }
455 theAngularDistribution[it]->SampleAndUpdat
456
457 if (theFinalStatePhotons[it] == nullptr) {
458 thePhotons = theGammas.GetDecayGammas(iL
459 eGamm -= theGammas.GetLevelEnergy(iLevel
460 }
461 }
462 else if (theEnergyAngData[it] != nullptr) /
463 {
464 theParticles = theEnergyAngData[it]->Sampl
465
466 // Adjust A and Z in the case of miss much
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)theParticl
473 auto ptr = theParticles->at(j);
474 G4int barnum = ptr->GetDefinition()->GetBary
475 if (barnum > maxA) {
476 maxA = barnum;
477 jAtMaxA = j;
478 }
479 sumA += barnum;
480 sumZ += G4lrint(ptr->GetDefinition()->
481 }
482 G4int dA = theBaseA + hadProjectile->Get
483 G4int dZ = theBaseZ +
484 G4lrint(hadProjectile->GetDefinition()->GetP
485 if (dA < 0 || dZ < 0) {
486 G4int newA = theParticles->at(jAtMaxA)
487 G4int newZ =
488 G4lrint(theParticles->at(jAtMaxA)->GetDefi
489 G4ParticleDefinition* pd = ionTable->G
490 theParticles->at(jAtMaxA)->SetDefiniti
491 }
492 }
493 }
494 else {
495 // @@@ what to do, if we have photon data,
496 nothingWasKnownOnHadron = 1;
497 }
498
499 if (theFinalStatePhotons[it] != nullptr) {
500 // the photon distributions are in the Nuc
501 // TK residual rest frame
502 G4ReactionProduct boosted_tmp;
503 boosted_tmp.Lorentz(incidReactionProduct,
504 G4double anEnergy = boosted_tmp.GetKinetic
505 thePhotons = theFinalStatePhotons[it]->Get
506 G4double aBaseEnergy = theFinalStatePhoton
507 G4double testEnergy = 0;
508 if (thePhotons != nullptr && !thePhotons->
509 aBaseEnergy -= (*thePhotons)[0]->GetTota
510 }
511 if (theFinalStatePhotons[it]->NeedsCascade
512 while (aBaseEnergy > 0.01 * CLHEP::keV)
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] != nullp
520 testEnergy = theFinalStatePhotons[
521 }
522 else {
523 testEnergy = 0;
524 }
525 G4double deltaE = std::abs(testEnerg
526 if (deltaE < 0.1 * CLHEP::keV) {
527 G4ReactionProductVector* theNext =
528 if (thePhotons != nullptr) thePhot
529 aBaseEnergy = testEnergy - theNext
530 delete theNext;
531 foundMatchingLevel = true;
532 break; // ===>
533 }
534 if (theFinalStatePhotons[j] != nullp
535 closest = j;
536 deltaEold = deltaE;
537 }
538 } // <=== the break goes here.
539 if (!foundMatchingLevel) {
540 G4ReactionProductVector* theNext = t
541 if (thePhotons != nullptr) thePhoton
542 aBaseEnergy = aBaseEnergy - theNext-
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 isotropi
557 // Next 12 lines are Emilio's replacement
558 // G4double QM=(incidReactionProduct.GetMa
559 // G4double eExcitation = QM-QI[it];
560 // G4double eExcitation = QI[0] - QI[it];
561 // if(eExcitation<20*CLHEP::keV){eExcitati
562
563 G4double eExcitation = std::max(0., QI[0]
564
565 two_body_reaction(&incidReactionProduct, &
566 if (thePhotons == nullptr && eExcitation >
567 for (iLevel = imaxEx; iLevel >= 0; --iLe
568 if (theGammas.GetLevelEnergy(iLevel) <
569 }
570 thePhotons = theGammas.GetDecayGammas(iL
571 }
572 }
573
574 // fill the result
575 // Beware - the recoil is not necessarily in
576 // Can be calculated from momentum conservat
577 // The idea is that the particles ar emitted
578 // recoil is on the residual; assumption is
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 < theParticl
589 aDef = (*theParticles)[ii0]->GetDefiniti
590 totalBaryonNumber += aDef->GetBaryonNumb
591 totalCharge += G4lrint(aDef->GetPDGCharg
592 totalMomentum += (*theParticles)[ii0]->G
593 }
594 if (totalBaryonNumber
595 != theBaseA + hadProjectile->GetDefin
596 {
597 needsSeparateRecoil = true;
598 residualA = theBaseA + hadProjectile->Ge
599 - totalBaryonNumber;
600 residualZ = theBaseZ +
601 G4lrint((hadProjectile->GetDefinition()->Get
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.GetDefinitio
615 theSec->SetMomentum(aHadron.GetMomentum())
616 theResult.Get()->AddSecondary(theSec, secI
617 #ifdef G4VERBOSE
618 if (fManager->GetDEBUG())
619 G4cout << " G4ParticleHPInelasticCompFS:
620 << theSec->GetParticleDefinition(
621 << " E= " << theSec->GetKineticEn
622 << theResult.Get()->GetNumberOfSe
623 #endif
624
625 aHadron.Lorentz(aHadron, theTarget);
626 G4ReactionProduct theResidual;
627 theResidual.SetDefinition(ionTable->GetIon
628 theResidual.SetKineticEnergy(aHadron.GetKi
629 / theResidual
630
631 // 080612TK contribution from Benoit Pirar
632 // theResidual.SetMomentum(-1.*aHadron.Get
633 G4ThreeVector incidentNeutronMomentum = in
634 theResidual.SetMomentum(incidentNeutronMom
635
636 theResidual.Lorentz(theResidual, -1. * the
637 G4ThreeVector totalPhotonMomentum(0, 0, 0)
638 if (thePhotons != nullptr) {
639 for (std::size_t i = 0; i < nPhotons; ++
640 totalPhotonMomentum += (*thePhotons)[i
641 }
642 }
643 theSec = new G4DynamicParticle;
644 theSec->SetDefinition(theResidual.GetDefin
645 theSec->SetMomentum(theResidual.GetMomentu
646 theResult.Get()->AddSecondary(theSec, secI
647 #ifdef G4VERBOSE
648 if (fManager->GetDEBUG())
649 G4cout << this << " G4ParticleHPInelasti
650 << theSec->GetParticleDefinition(
651 << " E= " << theSec->GetKineticEn
652 << theResult.Get()->GetNumberOfSe
653 #endif
654 }
655 else {
656 for (std::size_t i0 = 0; i0 < theParticles
657 theSec = new G4DynamicParticle;
658 theSec->SetDefinition((*theParticles)[i0
659 theSec->SetMomentum((*theParticles)[i0]-
660 theResult.Get()->AddSecondary(theSec, se
661 #ifdef G4VERBOSE
662 if (fManager->GetDEBUG())
663 G4cout << " G4ParticleHPInelasticCompF
664 << theSec->GetParticleDefinitio
665 << " E= " << theSec->GetKinetic
666 << theResult.Get()->GetNumberOf
667 #endif
668 delete (*theParticles)[i0];
669 }
670 delete theParticles;
671 if (needsSeparateRecoil && residualZ != 0)
672 G4ReactionProduct theResidual;
673 theResidual.SetDefinition(ionTable->GetI
674 G4double resiualKineticEnergy = theResid
675 resiualKineticEnergy += totalMomentum *
676 resiualKineticEnergy = std::sqrt(resiual
677 theResidual.SetKineticEnergy(resiualKine
678
679 // 080612TK contribution from Benoit Pir
680 // theResidual.SetMomentum(-1.*totalMome
681 // G4ThreeVector incidentNeutronMomentum
682 // theResidual.SetMomentum(incidentNeutr
683 // 080717 TK Comment still do NOT includ
684 theResidual.SetMomentum(incidReactionPro
685 - totalMomentum)
686
687 theSec = new G4DynamicParticle;
688 theSec->SetDefinition(theResidual.GetDef
689 theSec->SetMomentum(theResidual.GetMomen
690 theResult.Get()->AddSecondary(theSec, se
691 #ifdef G4VERBOSE
692 if (fManager->GetDEBUG())
693 G4cout << " G4ParticleHPInelasticCompF
694 << theSec->GetParticleDefinitio
695 << " E= " << theSec->GetKinetic
696 << theResult.Get()->GetNumberOf
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@u
704 // 2009 theSec->SetDefinition(G4Gamma::G
705 theSec->SetDefinition((*thePhotons)[i]->
706 // But never cause real effect at least
707 theSec->SetMomentum((*thePhotons)[i]->Ge
708 theResult.Get()->AddSecondary(theSec, se
709 #ifdef G4VERBOSE
710 if (fManager->GetDEBUG())
711 G4cout << " G4ParticleHPInelasticCompF
712 << theSec->GetParticleDefinitio
713 << " E= " << theSec->GetKinetic
714 << theResult.Get()->GetNumberOf
715 #endif
716
717 delete thePhotons->operator[](i);
718 }
719 // some garbage collection
720 delete thePhotons;
721 }
722
723 G4ParticleDefinition* targ_pd = ionTable->Ge
724 G4LorentzVector targ_4p_lab(
725 theTarget.GetMomentum(),
726 std::sqrt(targ_pd->GetPDGMass() * targ_pd-
727 G4LorentzVector proj_4p_lab = theTrack.Get4M
728 G4LorentzVector init_4p_lab = proj_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). Isotro
736 // proj: projectile in target-rest-frame (inp
737 // targ: target in target-rest-frame (input)
738 // product: secondary particle in target-rest
739 // resExcitationEnergy: excitation energy of
740 //
741 void G4ParticleHPInelasticCompFS::two_body_rea
742
743
744
745 {
746 // CMS system:
747 G4ReactionProduct theCMS = *proj + *targ;
748
749 // Residual definition:
750 G4int resZ = G4lrint((proj->GetDefinition()-
751 - product->GetDefinition()->GetPDGCharge
752 G4int resA = proj->GetDefinition()->GetBaryo
753 - product->GetDefinition()->Get
754 G4ReactionProduct theResidual;
755 theResidual.SetDefinition(ionTable->GetIon(r
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
764 + theCMSproj.GetTo
765 + std::sqrt(theCMStarg.GetMa
766 + theCMStarg.Get
767 G4double prodmass = product->GetMass();
768 G4double resmass = theResidual.GetMass() + r
769 G4double fmomsquared = (totE * totE - (prodm
770 (totE * totE - (prodmass + resmass) * (pro
771 G4double fmom = (fmomsquared > 0) ? std::sqr
772
773 // random (isotropic direction):
774 product->SetMomentum(fmom * G4RandomDirectio
775 product->SetTotalEnergy(std::sqrt(prodmass *
776 // Back to the LAB system:
777 product->Lorentz(*product, -1. * theCMS);
778 }
779
780 G4bool G4ParticleHPInelasticCompFS::use_nresp7
781
782
783
784 {
785 if (aDefinition == G4Neutron::Definition())
786 {
787 // LR: flag LR in ENDF. It indicates wheth
788 // it: exit channel (index of the carbon e
789
790 // Added by A. R. Garcia (CIEMAT) to inclu
791
792 if (LR[itt] > 0) // If there is breakup of
793 // MT=52-91 (it=MT-50)).
794 {
795 // Defining carbon as the target in the
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 mechanism
805 if (itt == 41) {
806 // QI=QM=-7.275 MeV for C-0(N,N')C-C(3
807 // This is not the value of the QI of
808 // to the model. So we don't take it.
809 // we have calculated: QI=(mn+m12C)-(m
810 nresp71_model.ApplyMechanismI_NBeA2A(b
811 // N+C --> A[0]+9BE* | 9BE* --> N[1]+8
812 }
813 else {
814 nresp71_model.ApplyMechanismII_ACN2A(b
815 // N+C --> N'[0]+C* | C* --> A[1]+8BE
816 }
817
818 // Returning to the reference frame wher
819 for (auto& theProd : theProds) {
820 theProd.Lorentz(theProd, -1. * theTarg
821 theResult.Get()->AddSecondary(
822 new G4DynamicParticle(theProd.GetDef
823 }
824
825 // Killing the primary neutron.
826 theResult.Get()->SetStatusChange(stopAnd
827
828 return true;
829 }
830 }
831 else if (aDefinition == G4Alpha::Definition(
832 {
833 // Added by A. R. Garcia (CIEMAT) to inclu
834
835 if (LR[itt] == 0) // If Z=6, an alpha part
836 // residual nucleus LR(f
837 {
838 // Defining carbon as the target in the
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,
848 // N+C --> A[0]+9BE[1].
849
850 for (auto& theProd : theProds) {
851 // Returning to the system of referenc
852 theProd.Lorentz(theProd, -1. * theTarg
853 theResult.Get()->AddSecondary(
854 new G4DynamicParticle(theProd.GetDef
855 }
856
857 // Killing the primary neutron.
858 theResult.Get()->SetStatusChange(stopAnd
859 return true;
860 }
861 G4Exception("G4ParticleHPInelasticCompFS::
862 FatalException, "Alpha product
863 }
864 return false;
865 }
866