Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4ParticleHPEnAngCorrelation.cc
Go to the documentation of this file.
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// 100413 Fix bug in incidence energy by T. Koi
31//
32// P. Arce, June-2014 Conversion neutron_hp to particle_hp
33//
34// June-2019 - E. Mendoza --> Part of the code trying to preserve the baryonic number has been
35// deleted. One has to assume that it is not preserved when using ENDF-6 data and it caused
36// problems.
37//
38// V. Ivanchenko, May-2023 Basic revision of particle HP classes
39//
40
42
43#include "G4IonTable.hh"
44#include "G4LorentzRotation.hh"
45#include "G4LorentzVector.hh"
46#include "G4RotationMatrix.hh"
48#include "G4Neutron.hh"
49
51{
52 theProjectile = (nullptr == proj) ? G4Neutron::Neutron() : proj;
53 toBeCached val;
54 fCache.Put( val );
55}
56
61
62void G4ParticleHPEnAngCorrelation::Init(std::istream& aDataFile)
63{
64 inCharge = true;
65 aDataFile >> targetMass >> frameFlag >> nProducts;
66 const std::size_t psize = nProducts > 0 ? nProducts : 1;
67 theProducts = new G4ParticleHPProduct[psize];
68 for (std::size_t i = 0; i < psize; ++i) {
69 theProducts[i].Init(aDataFile, theProjectile);
70 }
71}
72
74{
75 auto result = new G4ReactionProduct;
76
77 // do we have an appropriate distribution
78 if (nProducts != 1)
79 throw G4HadronicException(__FILE__, __LINE__, "More than one product in SampleOne");
80
81 // get the result
82 G4ReactionProductVector* temp = nullptr;
83 G4int i = 0;
84
85 G4int icounter = 0;
86 G4int icounter_max = 1024;
87 while (temp == nullptr)
88 {
89 ++icounter;
90 if (icounter > icounter_max)
91 {
92 G4cout << "Loop-counter exceeded the threshold value at "
93 << __LINE__ << "th line of "
94 << __FILE__ << "." << G4endl;
95 break;
96 }
97 temp = theProducts[++i].Sample(anEnergy, 1);
98 }
99 if (nullptr != temp)
100 {
101 if (temp->size() == 1)
102 {
103 result = (*temp)[0];
104 }
105 else
106 {
107 for (auto const& ptr : *temp) { delete ptr; }
108 throw G4HadronicException(__FILE__, __LINE__,
109 "SampleOne: Yield not correct");
110 }
111 }
112 delete temp;
113 return result;
114}
115
117{
118 auto result = new G4ReactionProductVector;
119 G4int i;
121 G4ReactionProduct theCMS;
123 /*
124 G4cout << "G4ParticleHPEnAngCorrelation::Sample for E=" << anEnergy
125 << " flag=" << frameFlag << " nProducts=" << nProducts <<G4endl;
126 */
127 if (frameFlag == 2 || frameFlag == 3) // Added for particle HP
128 {
129 G4ThreeVector the3IncidentPart = fCache.Get().theProjectileRP->GetMomentum();
130 // theProjectileRP has value in LAB
131 G4double nEnergy = fCache.Get().theProjectileRP->GetTotalEnergy();
132 G4ThreeVector the3Target = fCache.Get().theTarget->GetMomentum();
133 // theTarget has value in LAB
134 G4double tEnergy = fCache.Get().theTarget->GetTotalEnergy();
135 G4double totE = nEnergy + tEnergy;
136 G4ThreeVector the3CMS = the3Target + the3IncidentPart;
137 theCMS.SetMomentum(the3CMS);
138 G4double cmsMom = std::sqrt(the3CMS * the3CMS);
139 G4double sqrts = std::sqrt((totE - cmsMom) * (totE + cmsMom));
140 theCMS.SetMass(sqrts);
141 theCMS.SetTotalEnergy(totE);
142 G4ReactionProduct aIncidentPart;
143 aIncidentPart.Lorentz(*fCache.Get().theProjectileRP, theCMS);
144 // TKDB 100413
145 // ENDF-6 Formats Manual ENDF-102
146 // CHAPTER 6. FILE 6: PRODUCT ENERGY-ANGLE DISTRIBUTIONS
147 // LCT Reference system for secondary energy and angle
148 // incident energy is always given in the LAB system
149 anEnergy = fCache.Get().theProjectileRP->GetKineticEnergy();
150 // should be same argumment of "anEnergy"
151
152 G4LorentzVector Ptmp(aIncidentPart.GetMomentum(),
153 aIncidentPart.GetTotalEnergy());
154
155 toZ.rotateZ(-1 * Ptmp.phi());
156 toZ.rotateY(-1 * Ptmp.theta());
157 }
158 fCache.Get().theTotalMeanEnergy = 0;
159 G4LorentzRotation toLab(toZ.inverse());
160 // toLab only change axis NOT to LAB system
161
162 for (i = 0; i < nProducts; ++i) {
163 G4int nPart = theProducts[i].GetMultiplicity(anEnergy);
164 it = theProducts[i].Sample(anEnergy, nPart);
165 G4double aMeanEnergy = theProducts[i].MeanEnergyOfThisInteraction();
166
167 // 151120 TK Modified for solving reproducibility problem
168 // This change may have side effect.
169 if (aMeanEnergy >= 0) {
170 fCache.Get().theTotalMeanEnergy += aMeanEnergy;
171 }
172 else {
173 fCache.Get().theTotalMeanEnergy = anEnergy/nProducts + theProducts[i].GetQValue();
174 }
175 if (it != nullptr) {
176 for (unsigned int ii = 0; ii < it->size(); ++ii) {
177 G4LorentzVector pTmp1(it->operator[](ii)->GetMomentum(),
178 it->operator[](ii)->GetTotalEnergy());
179 pTmp1 = toLab * pTmp1;
180 it->operator[](ii)->SetMomentum(pTmp1.vect());
181 it->operator[](ii)->SetTotalEnergy(pTmp1.e());
182
183 if (frameFlag == 1) // target rest //TK 100413 should be LAB?
184 {
185 it->operator[](ii)->Lorentz(
186 *(it->operator[](ii)),
187 -1. * (*fCache.Get().theTarget)); // TK 100413 Is this really need?
188 }
189 else if (frameFlag == 2) // CMS
190 {
191#ifdef G4PHPDEBUG
192 if (G4ParticleHPManager::GetInstance()->GetDEBUG())
193 G4cout << "G4ParticleHPEnAngCorrelation: before Lorentz boost "
194 << it->at(ii)->GetKineticEnergy() << " "
195 << it->at(ii)->GetMomentum() << G4endl;
196#endif
197 it->operator[](ii)->Lorentz(*(it->operator[](ii)), -1. * theCMS);
198#ifdef G4PHPDEBUG
199 if (G4ParticleHPManager::GetInstance()->GetDEBUG())
200 G4cout << "G4ParticleHPEnAngCorrelation: after Lorentz boost "
201 << it->at(ii)->GetKineticEnergy() << " "
202 << it->at(ii)->GetMomentum() << G4endl;
203#endif
204 }
205 // TK120515 migrate frameFlag (MF6 LCT) = 3
206 else if (frameFlag == 3) // CMS A<=4 other LAB
207 {
208 if (theProducts[i].GetMassCode() > 2004.5) // Alpha AWP 3.96713
209 {
210 // LAB
211 it->operator[](ii)->Lorentz(
212 *(it->operator[](ii)),
213 -1. * (*fCache.Get().theTarget)); // TK 100413 Is this really need?
214#ifdef G4PHPDEBUG
215 if (G4ParticleHPManager::GetInstance()->GetDEBUG())
216 G4cout << "G4ParticleHPEnAngCorrelation: after Lorentz boost "
217 << it->at(ii)->GetKineticEnergy() << " "
218 << it->at(ii)->GetMomentum() << G4endl;
219#endif
220 }
221 else {
222 // CMS
223 it->operator[](ii)->Lorentz(*(it->operator[](ii)), -1. * theCMS);
224#ifdef G4PHPDEBUG
225 if (G4ParticleHPManager::GetInstance()->GetDEBUG())
226 G4cout << "G4ParticleHPEnAngCorrelation: after Lorentz boost "
227 << it->at(ii)->GetKineticEnergy() << " "
228 << it->at(ii)->GetMomentum() << G4endl;
229#endif
230 }
231 }
232 else
233 {
234 throw G4HadronicException(__FILE__, __LINE__,
235 "Sample: The frame of the final state is not specified");
236 }
237 result->push_back(it->operator[](ii));
238 }
239 delete it;
240 }
241 }
242 return result;
243}
std::vector< G4ReactionProduct * > G4ReactionProductVector
double G4double
Definition G4Types.hh:83
int G4int
Definition G4Types.hh:85
#define G4endl
Definition G4ios.hh:67
G4GLOB_DLL std::ostream G4cout
HepLorentzRotation & rotateY(double delta)
HepLorentzRotation & rotateZ(double delta)
HepLorentzRotation inverse() const
double theta() const
Hep3Vector vect() const
value_type & Get() const
Definition G4Cache.hh:315
void Put(const value_type &val) const
Definition G4Cache.hh:321
static G4Neutron * Neutron()
Definition G4Neutron.cc:101
G4ReactionProduct * SampleOne(G4double anEnergy)
G4ParticleHPEnAngCorrelation(const G4ParticleDefinition *proj=nullptr)
G4ReactionProductVector * Sample(G4double anEnergy)
static G4ParticleHPManager * GetInstance()
G4ReactionProductVector * Sample(G4double anEnergy, G4int nParticles)
G4double MeanEnergyOfThisInteraction()
void Init(std::istream &aDataFile, const G4ParticleDefinition *projectile)
G4int GetMultiplicity(G4double anEnergy)
void SetMomentum(const G4double x, const G4double y, const G4double z)
void SetTotalEnergy(const G4double en)
G4double GetTotalEnergy() const
G4ThreeVector GetMomentum() const
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
void SetMass(const G4double mas)