| 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 // 080801 Give a warning message for irregular
31 // Introduce theNDLDataA,Z which has A
32 // 081024 G4NucleiPropertiesTable:: to G4Nucle
33 // 101111 Add Special treatment for Be9(n,2n)B
34 //
35 // P. Arce, June-2014 Conversion neutron_hp to
36 //
37 // June-2019 - E. Mendoza --> Added protection
38 // is not applied when data is in MF=6 format
39 // (add Q value info to G4ParticleHPNBodyPhase
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::G4ParticleHPInela
53 {
54 hasXsec = true;
55 theXsection = new G4ParticleHPVector;
56 }
57
58 G4ParticleHPInelasticBaseFS::~G4ParticleHPInel
59 {
60 delete theXsection;
61 delete theEnergyDistribution;
62 delete theFinalStatePhotons;
63 delete theEnergyAngData;
64 delete theAngularDistribution;
65 }
66
67 void G4ParticleHPInelasticBaseFS::InitGammas(G
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 thi
77 std::ifstream theGammaData(aName, std::ios::
78
79 theNuclearMassDifference = G4NucleiPropertie
80 - G4NucleiProperties::GetBindingEnergy(the
81 theGammas.Init(theGammaData);
82 }
83
84 void G4ParticleHPInelasticBaseFS::Init(G4doubl
85 const G
86 const G
87 {
88 gammaPath = fManager->GetNeutronHPPath() + "
89 G4String tString = dirName;
90 SetA_Z(A, Z, M);
91 G4bool dbool = true;
92 G4ParticleHPDataUsed aFile =
93 theNames.GetName(theBaseA, theBaseZ, M, tS
94 SetAZMs(aFile);
95 G4String filename = aFile.GetName();
96 #ifdef G4VERBOSE
97 if (fManager->GetDEBUG())
98 G4cout << " G4ParticleHPInelasticBaseFS::I
99 #endif
100
101 if (!dbool || (theBaseZ <= 2 && (theNDLDataZ
102 {
103 #ifdef G4PHPDEBUG
104 if (fManager->GetDEBUG())
105 G4cout << "Skipped = " << filename << "
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 isoto
122 }
123 // here we go
124 G4int infoType, dataType, dummy = INT_MAX;
125 hasFSData = false;
126 while (theData >> infoType) // Loop checkin
127 {
128 theData >> dataType;
129
130 if (dummy == INT_MAX) theData >> Qvalue >>
131 Qvalue *= CLHEP::eV;
132 // In G4NDL4.5 this value is the MT number
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:
139 }
140 else if (dataType == 4) {
141 theAngularDistribution = new G4ParticleH
142 theAngularDistribution->Init(theData);
143 hasFSData = true;
144 }
145 else if (dataType == 5) {
146 theEnergyDistribution = new G4ParticleHP
147 theEnergyDistribution->Init(theData);
148 hasFSData = true;
149 }
150 else if (dataType == 6) {
151 theEnergyAngData = new G4ParticleHPEnAng
152 theEnergyAngData->Init(theData);
153 hasFSData = true;
154 }
155 else if (dataType == 12) {
156 theFinalStatePhotons = new G4ParticleHPP
157 theFinalStatePhotons->InitMean(theData);
158 hasFSData = true;
159 }
160 else if (dataType == 13) {
161 theFinalStatePhotons = new G4ParticleHPP
162 theFinalStatePhotons->InitPartials(theDa
163 hasFSData = true;
164 }
165 else if (dataType == 14) {
166 theFinalStatePhotons->InitAngular(theDat
167 hasFSData = true;
168 }
169 else if (dataType == 15) {
170 theFinalStatePhotons->InitEnergies(theDa
171 hasFSData = true;
172 }
173 else {
174 G4ExceptionDescription ed;
175 ed << "Z=" << theBaseZ << " A=" << theBa
176 << " projectile: " << theProjectile->GetPar
177 G4Exception("G4ParticleHPInelasticBaseFS
178 ed, "Data-type unknown");
179 }
180 }
181 }
182
183 void G4ParticleHPInelasticBaseFS::BaseApply(co
184 G4
185 G4
186 {
187 // G4cout << "G4ParticleHPInelasticBaseFS::B
188 // prepare neutron
189 if (theResult.Get() == nullptr) theResult.Pu
190 theResult.Get()->Clear();
191 G4double eKinetic = theTrack.GetKineticEnerg
192 G4ReactionProduct incidReactionProduct(theTr
193 incidReactionProduct.SetMomentum(theTrack.Ge
194 incidReactionProduct.SetKineticEnergy(eKinet
195
196 // prepare target
197 G4double targetMass = G4NucleiProperties::Ge
198 /CLHEP::neutron_mass_c2;
199 /*
200 G4cout << " Target Z=" << theBaseZ << " A="
201 << " projectile: " << theTrack.GetDefinitio
202 << " Ekin(MeV)=" << theTrack.GetKineticEner
203 << " Mtar/Mn=" << targetMass
204 << " hasFS=" << HasFSData() << G4endl;
205 */
206 // give priority to ENDF vales for target ma
207 if (theEnergyAngData != nullptr) {
208 G4double x = theEnergyAngData->GetTargetMa
209 if(0.0 < x) { targetMass = x/CLHEP::neutro
210 }
211 if (theAngularDistribution != nullptr) {
212 G4double x = theAngularDistribution->GetTa
213 if(0.0 < x) { targetMass = x/CLHEP::neutr
214 }
215
216 // 110512 TKDB ENDF-VII.0 21Sc45 has trouble
217 // recorded. TKDB targetMass = 0; ENDF-VII.0
218
219 G4Nucleus aNucleus;
220 G4ReactionProduct theTarget;
221
222 G4ThreeVector neuVelo =
223 incidReactionProduct.GetMomentum() / CLHEP
224
225 theTarget =
226 aNucleus.GetBiasedThermalNucleus(targetMas
227 theTrack.GetMaterial()->GetTemper
228 theTarget.SetDefinition(ionTable->GetIon(the
229
230 // prepare energy in target rest frame
231 G4ReactionProduct boosted;
232 boosted.Lorentz(incidReactionProduct, theTar
233 eKinetic = boosted.GetKineticEnergy();
234
235 G4double orgMomentum = boosted.GetMomentum()
236
237 // Take N-body phase-space distribution, if
238 if (!HasFSData()) // adding the residual is
239 {
240 G4ParticleHPNBodyPhaseSpace thePhaseSpaceD
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 calculate
250
251 // Calculate residual:
252 G4int ResidualA = theBaseA;
253 G4int ResidualZ = theBaseZ;
254 for (ii = 0; ii < nDef; ++ii) {
255 ResidualZ -= theDefs[ii]->GetAtomicNum
256 ResidualA -= theDefs[ii]->GetBaryonNum
257 }
258
259 if (ResidualA > 0 && ResidualZ > 0) {
260 G4ParticleDefinition* resid = ionTable
261 Qvalue =
262 incidReactionProduct.GetMass() + the
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,
273 thePhaseSpaceDistribution.SetProjectileRP(
274 thePhaseSpaceDistribution.SetTarget(&theTa
275 thePhaseSpaceDistribution.SetQValue(Qvalue
276
277 for (ii = 0; ii < nDef; ++ii) {
278 G4double massCode = 1000. * std::abs(the
279 massCode += theDefs[ii]->GetBaryonNumber
280 G4double dummy = 0;
281 G4ReactionProduct* aSec = thePhaseSpaceD
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, sec
288 #ifdef G4VERBOSE
289 if (fManager->GetDEBUG())
290 G4cout << this << " G4ParticleHPInelas
291 << aPart->GetParticleDefinition
292 << " E= " << aPart->GetKineticE
293 << theResult.Get()->GetNumberOf
294 #endif
295 }
296
297 theResult.Get()->SetStatusChange(stopAndKi
298 // Final momentum check should be done bef
299 const G4ParticleDefinition* targ_pd = ionT
300 G4double mass = targ_pd->GetPDGMass();
301 G4LorentzVector targ_4p_lab(theTarget.GetM
302 std::sqrt(mass*mass + theTarget.GetMo
303 G4LorentzVector proj_4p_lab = theTrack.Get
304 G4LorentzVector init_4p_lab = proj_4p_lab
305 adjust_final_state(init_4p_lab);
306 return;
307 }
308
309 // set target and neutron in the relevant ex
310 if (theAngularDistribution != nullptr) {
311 theAngularDistribution->SetTarget(theTarge
312 theAngularDistribution->SetProjectileRP(in
313 }
314 else if (theEnergyAngData != nullptr) {
315 theEnergyAngData->SetTarget(theTarget);
316 theEnergyAngData->SetProjectileRP(incidRea
317 }
318
319 G4ReactionProductVector* tmpHadrons = nullpt
320 G4int dummy;
321 std::size_t i;
322
323 if (theEnergyAngData != nullptr) {
324 tmpHadrons = theEnergyAngData->Sample(eKin
325 auto hproj = theTrack.GetDefinition();
326 G4int projA = hproj->GetBaryonNumber();
327 G4int projZ = G4lrint(hproj->GetPDGCharge(
328
329 if (!fManager->GetDoNotAdjustFinalState())
330 // Adjust A and Z in the case of miss mu
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)tmpHadro
337 auto had = ((*tmpHadrons)[j])->GetDefiniti
338 if (had->GetBaryonNumber() > maxA) {
339 maxA = had->GetBaryonNumber();
340 jAtMaxA = j;
341 }
342 sumA += had->GetBaryonNumber();
343 sumZ += G4lrint(had->GetPDGCharge()/
344 }
345 G4int dA = theBaseA + projA - sumA;
346 G4int dZ = theBaseZ + projZ - sumZ;
347 if (dA < 0 || dZ < 0) {
348 auto p = ((*tmpHadrons)[jAtMaxA])->G
349 G4int newA = p->GetBaryonNumber() +
350 G4int newZ = G4lrint(p->GetPDGCharge
351 if (newA > newZ && newZ > 0) {
352 G4ParticleDefinition* pd = ionTabl
353 ((*tmpHadrons)[jAtMaxA])->SetDefin
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::G
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(theEnergyD
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
384 }
385 else {
386 throw G4HadronicException(
387 __FILE__, __LINE__,
388 "No energy distribution to sample
389 }
390 theAngularDistribution->SampleAndUpdat
391 if (theEnergyDistribution == nullptr &
392 if (i0 == 0) {
393 G4double mass1 = theDefs[0]->GetPD
394 G4double mass2 = theDefs[1]->GetPD
395 G4double massn = theProjectile->Ge
396 G4int z1 = theBaseZ -
397 G4lrint((theDefs[0]->GetPDGCharge() +
398 G4int a1 = theBaseA - theDefs[0]->
399 G4double concreteMass = G4NucleiPr
400 G4double availableEnergy = eKineti
401 - concreteMass;
402 // available kinetic energy in CMS
403 G4double emin = availableEnergy +
404 - std::sqrt((mass1 + mass2) * (m
405 G4double p1 = std::sqrt(2. * mass2
406 bufferedDirection = p1 * aHadron->
407 #ifdef G4PHPDEBUG
408 if (fManager->GetDEBUG()) // @@@@
409 {
410 G4cout << "G4ParticleHPInelastic
411 << " " << a1 << " " << theBaseA << "
412 << " " << emin << G4endl;
413 }
414 #endif
415 }
416 else {
417 bufferedDirection = -bufferedDirec
418 }
419 // boost from cms to lab
420 aHadron->SetTotalEnergy(
421 std::sqrt(aHadron->GetMass() * aHa
422 aHadron->SetMomentum(bufferedDirecti
423 aHadron->Lorentz(*aHadron, -1. * (th
424 #ifdef G4PHPDEBUG
425 if (fManager->GetDEBUG()) {
426 G4cout << " G4ParticleHPInelasticB
427 << " P=" << aHadron->GetMomentum()
428 << " bufDir^2=" << bufferedDirection.ma
429 << G4endl;
430 }
431 #endif
432 }
433 tmpHadrons->push_back(aHadron);
434 #ifdef G4PHPDEBUG
435 if (fManager->GetDEBUG())
436 G4cout << " G4ParticleHPInelasticBas
437 << aHadron->GetDefinition()->
438 << " E= " << aHadron->GetKine
439 #endif
440 }
441 }
442 delete[] Done;
443 }
444 else {
445 throw G4HadronicException(__FILE__, __LINE
446 }
447
448 G4ReactionProductVector* thePhotons = nullpt
449 if (theFinalStatePhotons != nullptr) {
450 // the photon distributions are in the Nuc
451 G4ReactionProduct boosted_tmp;
452 boosted_tmp.Lorentz(incidReactionProduct,
453 G4double anEnergy = boosted_tmp.GetKinetic
454 thePhotons = theFinalStatePhotons->GetPhot
455 if (thePhotons != nullptr) {
456 for (i = 0; i < thePhotons->size(); ++i)
457 // back to lab
458 thePhotons->operator[](i)->Lorentz(*(t
459 }
460 }
461 }
462 else if (theEnergyAngData != nullptr) {
463 // PA130927: do not create photons to adju
464 G4bool bAdjustPhotons{true};
465 #ifdef PHP_AS_HP
466 bAdjustPhotons = true;
467 #else
468 if (fManager->GetDoNotAdjustFinalState())
469 #endif
470
471 if (bAdjustPhotons) {
472 G4double theGammaEnergy = theEnergyAngDa
473 G4double anEnergy = boosted.GetKineticEn
474 theGammaEnergy = anEnergy - theGammaEner
475 theGammaEnergy += theNuclearMassDifferen
476 G4double eBindProducts = 0;
477 G4double eBindN = 0;
478 G4double eBindP = 0;
479 G4double eBindD = G4NucleiProperties::Ge
480 G4double eBindT = G4NucleiProperties::Ge
481 G4double eBindHe3 = G4NucleiProperties::
482 G4double eBindA = G4NucleiProperties::Ge
483 G4int ia = 0;
484 for (auto const & had : *tmpHadrons) {
485 if (had->GetDefinition() == G4Neutron:
486 eBindProducts += eBindN;
487 }
488 else if (had->GetDefinition() == G4Pro
489 eBindProducts += eBindP;
490 }
491 else if (had->GetDefinition() == G4Deu
492 eBindProducts += eBindD;
493 }
494 else if (had->GetDefinition() == G4Tri
495 eBindProducts += eBindT;
496 }
497 else if (had->GetDefinition() == G4He3
498 eBindProducts += eBindHe3;
499 }
500 else if (had->GetDefinition() == G4Alp
501 eBindProducts += eBindA;
502 ia++;
503 }
504 }
505 theGammaEnergy += eBindProducts;
506
507 #ifdef G4PHPDEBUG
508 if (fManager->GetDEBUG())
509 G4cout << " G4ParticleHPInelasticBaseF
510 << " eBindProducts " << eBindPr
511 #endif
512
513 // Special treatment for Be9 + n -> 2n +
514 if (theBaseZ == 4 && theBaseA == 9) {
515 // This only valid for G4NDL3.13,,,
516 if (ia == 2 && std::abs(G4NucleiProper
517 G4NucleiProperties::GetBindingEnergy(9
518 theNuclearMassDifference) < CLHEP::keV
519 {
520 theGammaEnergy -= (2 * eBindA);
521 }
522 }
523
524 if (theGammaEnergy > 0.0) {
525 for (G4int iLevel = theGammas.GetNumberOfLev
526 G4double e = theGammas.GetLevelEnerg
527 if (e < theGammaEnergy) {
528 thePhotons = theGammas.GetDecayGam
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::B
548 << G4endl;
549 #endif
550
551 for (i = 0; i < nSecondaries - nPhotons; ++i
552 theSec = new G4DynamicParticle;
553 theSec->SetDefinition(tmpHadrons->operator
554 theSec->SetMomentum(tmpHadrons->operator[]
555 theResult.Get()->AddSecondary(theSec, secI
556 #ifdef G4PHPDEBUG
557 if (fManager->GetDEBUG())
558 G4cout << " G4ParticleHPInelasticBaseFS:
559 << theSec->GetParticleDefinition(
560 << " E=" << theSec->GetKineticEne
561 << theResult.Get()->GetNumberOfSe
562 #endif
563 delete (*tmpHadrons)[i];
564 }
565 #ifdef G4PHPDEBUG
566 if (fManager->GetDEBUG())
567 G4cout << " G4ParticleHPInelasticBaseFS::B
568 #endif
569 if (thePhotons != nullptr) {
570 for (i = 0; i < nPhotons; ++i) {
571 theSec = new G4DynamicParticle;
572 theSec->SetDefinition((*thePhotons)[i]->
573 theSec->SetMomentum((*thePhotons)[i]->Ge
574 theResult.Get()->AddSecondary(theSec, se
575 #ifdef G4PHPDEBUG
576 if (fManager->GetDEBUG())
577 G4cout << " G4ParticleHPInelasticBaseF
578 << theSec->GetParticleDefinitio
579 << " E=" << theSec->GetKineticE
580 << theResult.Get()->GetNumberOf
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->Ge
591 G4LorentzVector targ_4p_lab(
592 theTarget.GetMomentum(),
593 std::sqrt(targ_pd->GetPDGMass() * targ_pd-
594 G4LorentzVector proj_4p_lab = theTrack.Get4M
595 G4LorentzVector init_4p_lab = proj_4p_lab +
596
597 // if data in MF=6 format (no correlated par
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