| 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 // neutron_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 // 071031 bug fix T. Koi on behalf of A. Howar
32 // 081203 bug fix in Register method by T. Koi
33 //
34 // P. Arce, June-2014 Conversion neutron_hp to
35 //
36 // June-2019 - E. Mendoza --> Modification to
37 // data library if the G4NEUTRONHP_SKIP_MISS
38 // flag is defined. The missing XS are set t
39
40 #include "G4ParticleHPChannel.hh"
41
42 #include "G4HadTmpUtil.hh"
43 #include "G4ParticleHPElasticFS.hh"
44 #include "G4ParticleHPFinalState.hh"
45 #include "G4ParticleHPReactionWhiteBoard.hh"
46 #include "G4ParticleHPThermalBoost.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "globals.hh"
49
50 #include <cstdlib>
51
52 G4ParticleHPChannel::G4ParticleHPChannel(G4Par
53 {
54 fManager = G4ParticleHPManager::GetInstance(
55 if (fManager->GetUseWendtFissionModel()) {
56 wendtFissionGenerator = G4WendtFissionFrag
57 // Make sure both fission fragment models
58 fManager->SetProduceFissionFragments(false
59 }
60 theProjectile = (nullptr == p) ? G4Neutron::
61 theChannelData = new G4ParticleHPVector;
62 }
63
64 G4ParticleHPChannel::~G4ParticleHPChannel()
65 {
66 delete theChannelData;
67 // Following statement disabled to avoid SEG
68 // theBuffer is also deleted as "theChannelD
69 delete[] theIsotopeWiseData;
70 if (theFinalStates != nullptr) {
71 for (G4int i = 0; i < niso; i++) {
72 delete theFinalStates[i];
73 }
74 delete[] theFinalStates;
75 }
76 delete[] active;
77 }
78
79 G4double G4ParticleHPChannel::GetXsec(G4double
80 {
81 return std::max(0., theChannelData->GetXsec(
82 }
83
84 G4double G4ParticleHPChannel::GetWeightedXsec(
85 G4int isoNumber) const
86 {
87 return theIsotopeWiseData[isoNumber].GetXsec
88 }
89
90 G4double G4ParticleHPChannel::GetFSCrossSectio
91 G4int isoNumber) const
92 {
93 return theFinalStates[isoNumber]->GetXsec(en
94 }
95
96 void G4ParticleHPChannel::Init(G4Element* anEl
97 const G4String& dirName, const G4
98 {
99 theFSType = aFSType;
100 Init(anElement, dirName);
101 }
102
103 void G4ParticleHPChannel::Init(G4Element* anEl
104 {
105 theDir = dirName;
106 theElement = anElement;
107 }
108
109 G4bool G4ParticleHPChannel::Register(G4Particl
110 {
111 ++registerCount;
112 G4int Z = theElement->GetZasInt();
113
114 niso = (G4int)theElement->GetNumberOfIsotope
115 const std::size_t nsize = niso > 0 ? niso :
116
117 delete[] theIsotopeWiseData;
118 theIsotopeWiseData = new G4ParticleHPIsoData
119 delete[] active;
120 active = new G4bool[nsize];
121
122 delete[] theFinalStates;
123 theFinalStates = new G4ParticleHPFinalState*
124 delete theChannelData;
125 theChannelData = new G4ParticleHPVector;
126 for (G4int i = 0; i < niso; ++i) {
127 theFinalStates[i] = theFS->New();
128 theFinalStates[i]->SetProjectile(theProjec
129 }
130 if (niso != 0 && registerCount == 0) {
131 for (G4int i1 = 0; i1 < niso; ++i1) {
132 G4int A = theElement->GetIsotope(i1)->Ge
133 G4int M = theElement->GetIsotope(i1)->Ge
134 //G4cout <<" Init: normal case i=" << i1
135 // << " Z=" << Z << " A=" << A << G4
136 G4double frac = theElement->GetRelativeA
137 theFinalStates[i1]->SetA_Z(A, Z, M);
138 UpdateData(A, Z, M, i1, frac, theProject
139 }
140 }
141 G4bool result = HasDataInAnyFinalState();
142
143 // To avoid issuing hash by worker threads
144 if (result) theChannelData->Hash();
145
146 return result;
147 }
148
149 void G4ParticleHPChannel::UpdateData(G4int A,
150 G4double
151 G4Particl
152 {
153 // Initialze the G4FissionFragment generator
154 if (wendtFissionGenerator != nullptr) {
155 wendtFissionGenerator->InitializeANucleus(
156 }
157
158 theFinalStates[index]->Init(A, Z, M, theDir,
159 if (!theFinalStates[index]->HasAnyData()) re
160 // nothing there for exactly this isotope.
161
162 // the above has put the X-sec into the FS
163 theBuffer = nullptr;
164 if (theFinalStates[index]->HasXsec()) {
165 theBuffer = theFinalStates[index]->GetXsec
166 theBuffer->Times(abundance / 100.);
167 theIsotopeWiseData[index].FillChannelData(
168 }
169 else // get data from CrossSection director
170 {
171 const G4String& tString = "/CrossSection";
172 active[index] = theIsotopeWiseData[index].
173
174 if (active[index]) theBuffer = theIsotopeW
175 }
176 if (theBuffer != nullptr) Harmonise(theChann
177 }
178
179 void G4ParticleHPChannel::Harmonise(G4Particle
180 G4Particle
181 {
182 G4int s_tmp = 0, n = 0, m_tmp = 0;
183 auto theMerge = new G4ParticleHPVector;
184 G4ParticleHPVector* anActive = theStore;
185 G4ParticleHPVector* aPassive = theNew;
186 G4ParticleHPVector* tmp;
187 G4int a = s_tmp, p = n, t;
188 while (a < anActive->GetVectorLength() && p
189 // Loop checking, 11.05.2015, T. Koi
190 {
191 if (anActive->GetEnergy(a) <= aPassive->Ge
192 G4double xa = anActive->GetEnergy(a);
193 theMerge->SetData(m_tmp, xa, anActive->G
194 m_tmp++;
195 a++;
196 G4double xp = aPassive->GetEnergy(p);
197 if (std::abs(std::abs(xp - xa) / xa) < 0
198 ++p;
199 }
200 }
201 else {
202 tmp = anActive;
203 t = a;
204 anActive = aPassive;
205 a = p;
206 aPassive = tmp;
207 p = t;
208 }
209 }
210 while (a != anActive->GetVectorLength()) //
211 {
212 theMerge->SetData(m_tmp++, anActive->GetEn
213 ++a;
214 }
215 while (p != aPassive->GetVectorLength()) //
216 {
217 if (std::abs(theMerge->GetEnergy(std::max(
218 aPassive->GetEnergy(p)) / aPassive->GetEn
219 theMerge->SetData(m_tmp++, aPassive->Get
220 ++p;
221 }
222 delete theStore;
223 theStore = theMerge;
224 }
225
226 G4WendtFissionFragmentGenerator* G4ParticleHPC
227 if ( wendtFissionGenerator ) return wendtFis
228 else return nullptr;
229 }
230
231 G4HadFinalState*
232 G4ParticleHPChannel::ApplyYourself(const G4Had
233 G4int anIsotope, G4bool isElastic)
234 {
235 //G4cout << "G4ParticleHPChannel::ApplyYours
236 // << " ni=" << anIsotope << " isElastic="
237 if (anIsotope != -1 && anIsotope != -2) {
238 // Inelastic Case
239 //G4cout << "G4ParticleHPChannel Inelastic
240 //<< " Z= " << GetZ(anIsotope) << " A = "
241 fManager->GetReactionWhiteBoard()->SetTarg
242 fManager->GetReactionWhiteBoard()->SetTarg
243 return theFinalStates[anIsotope]->ApplyYou
244 }
245 G4double sum = 0;
246 G4int it = 0;
247 auto xsec = new G4double[niso];
248 G4ParticleHPThermalBoost aThermalE;
249 for (G4int i = 0; i < niso; i++) {
250 if (theFinalStates[i]->HasAnyData()) {
251 /*
252 G4cout << "FS: " << i << theTrack.GetDef
253 << " Z=" << theFinalStates[i]->GetZ()
254 << " A=" << theFinalStates[i]->GetN()
255 << G4endl;
256 */
257 xsec[i] = theIsotopeWiseData[i].GetXsec(
258 aThermalE.GetThermalEnergy(theTrack, t
259 theFinalSta
260 theTrack.Ge
261 sum += xsec[i];
262 }
263 else {
264 xsec[i] = 0;
265 }
266 }
267 if (sum == 0) {
268 it = G4lrint(niso * G4UniformRand());
269 }
270 else {
271 G4double random = G4UniformRand();
272 G4double running = 0;
273 for (G4int ix = 0; ix < niso; ix++) {
274 running += xsec[ix];
275 if (sum == 0 || random <= running / sum)
276 it = ix;
277 break;
278 }
279 }
280 if (it == niso) it--;
281 }
282 delete[] xsec;
283 G4HadFinalState* theFinalState = nullptr;
284 const auto A = (G4int)this->GetN(it);
285 const auto Z = (G4int)this->GetZ(it);
286 const auto M = (G4int)this->GetM(it);
287
288 //-2:Marker for Fission
289 if ((wendtFissionGenerator != nullptr) && an
290 theFinalState = wendtFissionGenerator->App
291 }
292
293 // Use the standard procedure if the G4Fissi
294 if (theFinalState == nullptr) {
295 G4int icounter = 0;
296 G4int icounter_max = 1024;
297 while (theFinalState == nullptr) // Loop
298 {
299 icounter++;
300 if (icounter > icounter_max) {
301 G4cout << "Loop-counter exceeded the t
302 << __LINE__ << "th line of " <<
303 break;
304 }
305 if (isElastic) {
306 // Register 0 K cross-section for DBRC
307 G4ParticleHPVector* xsforFS = theIsoto
308 // Only G4ParticleHPElasticFS has the
309 static_cast<G4ParticleHPElasticFS*>(th
310 }
311 theFinalState = theFinalStates[it]->Appl
312 }
313 }
314
315 // G4cout <<"THE IMPORTANT RETURN"<<G4endl;
316 // G4cout << "TK G4ParticleHPChannel Elastic
317 //<< " Z= " << this->GetZ(it) << " A = " <<
318 fManager->GetReactionWhiteBoard()->SetTargA(
319 fManager->GetReactionWhiteBoard()->SetTargZ(
320 fManager->GetReactionWhiteBoard()->SetTargM(
321
322 return theFinalState;
323 }
324
325 void G4ParticleHPChannel::DumpInfo() const
326 {
327 G4cout << " Element: " << theElement->GetNam
328 G4cout << " Directory name: " << theDir << G
329 G4cout << " FS name: " << theFSType << G4end
330 G4cout << " Number of Isotopes: " << niso <<
331 G4cout << " Have cross sections: " << G4endl
332 for (int i = 0; i < niso; i++) {
333 G4cout << theFinalStates[i]->HasXsec() <<
334 }
335 G4cout << G4endl;
336 if (theChannelData != nullptr) {
337 G4cout << " Cross Section (total for this
338 int np = theChannelData->GetVectorLength()
339 G4cout << np << G4endl;
340 for (int i = 0; i < np; i++) {
341 G4cout << theChannelData->GetEnergy(i) /
342 }
343 }
344 }
345