| 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 // 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
41 #include "G4ParticleHPEnAngCorrelation.hh"
42
43 #include "G4IonTable.hh"
44 #include "G4LorentzRotation.hh"
45 #include "G4LorentzVector.hh"
46 #include "G4RotationMatrix.hh"
47 #include "G4ParticleHPManager.hh"
48 #include "G4Neutron.hh"
49
50 G4ParticleHPEnAngCorrelation::G4ParticleHPEnAngCorrelation(const G4ParticleDefinition* proj)
51 {
52 theProjectile = (nullptr == proj) ? G4Neutron::Neutron() : proj;
53 toBeCached val;
54 fCache.Put( val );
55 }
56
57 G4ParticleHPEnAngCorrelation::~G4ParticleHPEnAngCorrelation()
58 {
59 delete[] theProducts;
60 }
61
62 void 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
73 G4ReactionProduct* G4ParticleHPEnAngCorrelation::SampleOne(G4double anEnergy)
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
116 G4ReactionProductVector* G4ParticleHPEnAngCorrelation::Sample(G4double anEnergy)
117 {
118 auto result = new G4ReactionProductVector;
119 G4int i;
120 G4ReactionProductVector* it;
121 G4ReactionProduct theCMS;
122 G4LorentzRotation toZ;
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 }
244