| 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 // 080808 Bug fix in serching mu bin and index
31 //
32 // P. Arce, June-2014 Conversion neutron_hp to
33 //
34 // E. Mendoza, Nov. 2020 - bug fix
35 //
36 #include "G4ParticleHPLabAngularEnergy.hh"
37
38 #include "G4Alpha.hh"
39 #include "G4Deuteron.hh"
40 #include "G4Electron.hh"
41 #include "G4Gamma.hh"
42 #include "G4He3.hh"
43 #include "G4Neutron.hh"
44 #include "G4PhysicalConstants.hh"
45 #include "G4Positron.hh"
46 #include "G4Proton.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4Triton.hh"
49 #include "Randomize.hh"
50
51 void G4ParticleHPLabAngularEnergy::Init(std::i
52 {
53 aDataFile >> nEnergies;
54 theManager.Init(aDataFile);
55 const std::size_t esize = nEnergies > 0 ? nE
56 theEnergies = new G4double[esize];
57 nCosTh = new G4int[esize];
58 theData = new G4ParticleHPVector*[esize];
59 theSecondManager = new G4InterpolationManage
60 for (G4int i = 0; i < nEnergies; ++i) {
61 aDataFile >> theEnergies[i];
62 theEnergies[i] *= eV;
63 aDataFile >> nCosTh[i];
64 theSecondManager[i].Init(aDataFile);
65 const std::size_t dsize = nCosTh[i] > 0 ?
66 theData[i] = new G4ParticleHPVector[dsize]
67 G4double label;
68 for (std::size_t ii = 0; ii < dsize; ++ii)
69 aDataFile >> label;
70 theData[i][ii].SetLabel(label);
71 theData[i][ii].Init(aDataFile, eV);
72 }
73 }
74 }
75
76 G4ReactionProduct* G4ParticleHPLabAngularEnerg
77
78 {
79 auto result = new G4ReactionProduct;
80 auto Z = static_cast<G4int>(massCode / 1000)
81 auto A = static_cast<G4int>(massCode - 1000
82
83 if (massCode == 0) {
84 result->SetDefinition(G4Gamma::Gamma());
85 }
86 else if (A == 0) {
87 result->SetDefinition(G4Electron::Electron
88 if (Z == 1) result->SetDefinition(G4Positr
89 }
90 else if (A == 1) {
91 result->SetDefinition(G4Neutron::Neutron()
92 if (Z == 1) result->SetDefinition(G4Proton
93 }
94 else if (A == 2) {
95 result->SetDefinition(G4Deuteron::Deuteron
96 }
97 else if (A == 3) {
98 result->SetDefinition(G4Triton::Triton());
99 if (Z == 2) result->SetDefinition(G4He3::H
100 }
101 else if (A == 4) {
102 result->SetDefinition(G4Alpha::Alpha());
103 if (Z != 2) throw G4HadronicException(__FI
104 }
105 else {
106 throw G4HadronicException(__FILE__, __LINE
107 "G4ParticleHPLab
108 }
109
110 // get theta, E
111 G4double cosTh, secEnergy;
112 G4int i, it(0);
113 // find the energy bin
114 for (i = 0; i < nEnergies; i++) {
115 it = i;
116 if (anEnergy < theEnergies[i]) break;
117 }
118
119 if (it == 0) // it marks the energy bin -->
120 // to high energies (??)
121 {
122 // Do not sample between incident neutron
123 //----------------------------------------
124 // CosTheta:
125 G4ParticleHPVector theThVec;
126 theThVec.SetInterpolationManager(theSecond
127 for (i = 0; i < nCosTh[it]; i++) {
128 theThVec.SetX(i, theData[it][i].GetLabel
129 theThVec.SetY(i, theData[it][i].GetInteg
130 }
131 cosTh = theThVec.Sample();
132 //----------------------------------------
133
134 //----------------------------------------
135 // Now the secondary energy:
136 G4double x, x1, x2, y1, y2, y, E;
137 G4int i1(0);
138 for (i = 1; i < nCosTh[it]; i++) {
139 i1 = i;
140 if (cosTh < theData[it][i].GetLabel()) b
141 }
142 // now get the prob at this energy for the
143 x = cosTh;
144 x1 = theData[it][i1 - 1].GetLabel();
145 x2 = theData[it][i1].GetLabel();
146 G4ParticleHPVector theBuff1a;
147 theBuff1a.SetInterpolationManager(theData[
148 for (i = 0; i < theData[it][i1 - 1].GetVec
149 E = theData[it][i1 - 1].GetX(i);
150 y1 = theData[it][i1 - 1].GetY(i);
151 y2 = theData[it][i1].GetY(E);
152 y = theInt.Interpolate(theSecondManager[
153 theBuff1a.SetData(i, E, y);
154 }
155 G4ParticleHPVector theBuff2a;
156 theBuff2a.SetInterpolationManager(theData[
157 for (i = 0; i < theData[it][i1].GetVectorL
158 E = theData[it][i1].GetX(i);
159 y1 = theData[it][i1 - 1].GetY(E);
160 y2 = theData[it][i1].GetY(i);
161 y = theInt.Interpolate(theSecondManager[
162 theBuff2a.SetData(i, E, y); // wrong E,
163 }
164 G4ParticleHPVector theStore;
165 theStore.Merge(&theBuff1a, &theBuff2a);
166 secEnergy = theStore.Sample();
167 currentMeanEnergy = theStore.GetMeanX();
168 //----------------------------------------
169 }
170 else // this is the small big else.
171 {
172 G4double x, x1, x2, y1, y2, y, tmp, E;
173 // integrate the prob for each costh, and
174 G4ParticleHPVector run1;
175 for (i = 0; i < nCosTh[it - 1]; i++) {
176 run1.SetX(i, theData[it - 1][i].GetLabel
177 run1.SetY(i, theData[it - 1][i].GetInteg
178 }
179 G4ParticleHPVector run2;
180 for (i = 0; i < nCosTh[it]; i++) {
181 run2.SetX(i, theData[it][i].GetLabel());
182 run2.SetY(i, theData[it][i].GetIntegral(
183 }
184 // get the distributions for the correct n
185 x = anEnergy;
186 x1 = theEnergies[it - 1];
187 x2 = theEnergies[it];
188 G4ParticleHPVector thBuff1; // to be inte
189 thBuff1.SetInterpolationManager(theSecondM
190 for (i = 0; i < run1.GetVectorLength(); i+
191 tmp = run1.GetX(i); // theta
192 y1 = run1.GetY(i); // integral
193 y2 = run2.GetY(tmp);
194 y = theInt.Interpolate(theManager.GetSch
195 thBuff1.SetData(i, tmp, y);
196 }
197 G4ParticleHPVector thBuff2;
198 thBuff2.SetInterpolationManager(theSecondM
199 for (i = 0; i < run2.GetVectorLength(); i+
200 tmp = run2.GetX(i); // theta
201 y1 = run1.GetY(tmp); // integral
202 y2 = run2.GetY(i);
203 y = theInt.Interpolate(theManager.GetSch
204 thBuff2.SetData(i, tmp, y);
205 }
206 G4ParticleHPVector theThVec;
207 theThVec.Merge(&thBuff1, &thBuff2); // ta
208 cosTh = theThVec.Sample();
209 /*
210 if(massCode>0.99 && massCode<1.01){theThVe
211 G4double random = (theThVec.GetY(theThVec.
212 -theThVec.GetY(0)) *G
213 G4cout<<" -- "<<theThVec.GetY(theThVec.Get
214 "<<random<<G4endl; G4int ith(0); for(i=1;i
215 {
216 ith = i;
217 if(random<theThVec.GetY(i)-theThVec.GetY
218 }
219 {
220 // calculate theta
221 G4double xx, xx1, xx2, yy1, yy2;
222 xx = random;
223 xx1 = theThVec.GetY(ith-1)-theThVec.GetY
224 xx2 = theThVec.GetY(ith)-theThVec.GetY(0
225 yy1 = theThVec.GetX(ith-1); // std::cos(
226 yy2 = theThVec.GetX(ith);
227 cosTh = theInt.Interpolate(theSecondMana
228 xx, xx1,xx2,y
229 }
230 */
231 G4int i1(0), i2(0);
232 // get the indixes of the vectors close to
233 // first it-1 !!!! i.e. low in energy
234 for (i = 1; i < nCosTh[it - 1]; i++) {
235 i1 = i;
236 if (cosTh < theData[it - 1][i].GetLabel(
237 }
238 // now get the prob at this energy for the
239 x = cosTh;
240 x1 = theData[it - 1][i1 - 1].GetLabel();
241 x2 = theData[it - 1][i1].GetLabel();
242 G4ParticleHPVector theBuff1a;
243 theBuff1a.SetInterpolationManager(theData[
244 for (i = 0; i < theData[it - 1][i1 - 1].Ge
245 E = theData[it - 1][i1 - 1].GetX(i);
246 y1 = theData[it - 1][i1 - 1].GetY(i);
247 y2 = theData[it - 1][i1].GetY(E);
248 y = theInt.Interpolate(theSecondManager[
249 theBuff1a.SetData(i, E, y); // wrong E,
250 }
251 G4ParticleHPVector theBuff2a;
252 theBuff2a.SetInterpolationManager(theData[
253 for (i = 0; i < theData[it - 1][i1].GetVec
254 E = theData[it - 1][i1].GetX(i);
255 y1 = theData[it - 1][i1 - 1].GetY(E);
256 y2 = theData[it - 1][i1].GetY(i);
257 y = theInt.Interpolate(theSecondManager[
258 theBuff2a.SetData(i, E, y); // wrong E,
259 }
260 G4ParticleHPVector theStore1;
261 theStore1.Merge(&theBuff1a, &theBuff2a);
262
263 // get the indixes of the vectors close to
264 // then it !!!! i.e. high in energy
265 for (i = 1; i < nCosTh[it]; i++) {
266 i2 = i;
267 if (cosTh < theData[it][i2].GetLabel())
268 } // sonderfaelle mit i1 oder i2 head on
269 x1 = theData[it][i2 - 1].GetLabel();
270 x2 = theData[it][i2].GetLabel();
271 G4ParticleHPVector theBuff1b;
272 theBuff1b.SetInterpolationManager(theData[
273 for (i = 0; i < theData[it][i2 - 1].GetVec
274 E = theData[it][i2 - 1].GetX(i);
275 y1 = theData[it][i2 - 1].GetY(i);
276 y2 = theData[it][i2].GetY(E);
277 y = theInt.Interpolate(theSecondManager[
278 theBuff1b.SetData(i, E, y); // wrong E,
279 }
280 G4ParticleHPVector theBuff2b;
281 theBuff2b.SetInterpolationManager(theData[
282 // 080808 i1 -> i2
283 // for(i=0;i<theData[it][i1].GetVectorLeng
284 for (i = 0; i < theData[it][i2].GetVectorL
285 // E = theData[it][i1].GetX(i);
286 // y1 = theData[it][i1-1].GetY(E);
287 // y2 = theData[it][i1].GetY(i);
288 E = theData[it][i2].GetX(i);
289 y1 = theData[it][i2 - 1].GetY(E);
290 y2 = theData[it][i2].GetY(i);
291 y = theInt.Interpolate(theSecondManager[
292 theBuff2b.SetData(i, E, y); // wrong E,
293 }
294 G4ParticleHPVector theStore2;
295 theStore2.Merge(&theBuff1b, &theBuff2b);
296 // now get to the right energy.
297
298 x = anEnergy;
299 x1 = theEnergies[it - 1];
300 x2 = theEnergies[it];
301 G4ParticleHPVector theOne1;
302 theOne1.SetInterpolationManager(theStore1.
303 for (i = 0; i < theStore1.GetVectorLength(
304 E = theStore1.GetX(i);
305 y1 = theStore1.GetY(i);
306 y2 = theStore2.GetY(E);
307 y = theInt.Interpolate(theManager.GetSch
308 theOne1.SetData(i, E, y); // both corre
309 }
310 G4ParticleHPVector theOne2;
311 theOne2.SetInterpolationManager(theStore2.
312 for (i = 0; i < theStore2.GetVectorLength(
313 E = theStore2.GetX(i);
314 y1 = theStore1.GetY(E);
315 y2 = theStore2.GetY(i);
316 y = theInt.Interpolate(theManager.GetSch
317 theOne2.SetData(i, E, y); // both corre
318 }
319 G4ParticleHPVector theOne;
320 theOne.Merge(&theOne1, &theOne2); // both
321
322 secEnergy = theOne.Sample();
323 currentMeanEnergy = theOne.GetMeanX();
324 }
325
326 // now do random direction in phi, and fill
327
328 result->SetKineticEnergy(secEnergy);
329
330 G4double phi = twopi * G4UniformRand();
331 G4double theta = std::acos(cosTh);
332 G4double sinth = std::sin(theta);
333 G4double mtot = result->GetTotalMomentum();
334 G4ThreeVector tempVector(mtot * sinth * std:
335 mtot * std::cos(the
336 result->SetMomentum(tempVector);
337
338 return result;
339 }
340