| 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 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron transport model.
29 //
30 // 080808 Bug fix in serching mu bin and index for theBuff2b by T. Koi
31 //
32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
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::istream& aDataFile)
52 {
53 aDataFile >> nEnergies;
54 theManager.Init(aDataFile);
55 const std::size_t esize = nEnergies > 0 ? nEnergies : 1;
56 theEnergies = new G4double[esize];
57 nCosTh = new G4int[esize];
58 theData = new G4ParticleHPVector*[esize];
59 theSecondManager = new G4InterpolationManager[esize];
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 ? nCosTh[i] : 1;
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* G4ParticleHPLabAngularEnergy::Sample(G4double anEnergy, G4double massCode,
77 G4double)
78 {
79 auto result = new G4ReactionProduct;
80 auto Z = static_cast<G4int>(massCode / 1000);
81 auto A = static_cast<G4int>(massCode - 1000 * Z);
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(G4Positron::Positron());
89 }
90 else if (A == 1) {
91 result->SetDefinition(G4Neutron::Neutron());
92 if (Z == 1) result->SetDefinition(G4Proton::Proton());
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::He3());
100 }
101 else if (A == 4) {
102 result->SetDefinition(G4Alpha::Alpha());
103 if (Z != 2) throw G4HadronicException(__FILE__, __LINE__, "Unknown ion case 1");
104 }
105 else {
106 throw G4HadronicException(__FILE__, __LINE__,
107 "G4ParticleHPLabAngularEnergy: Unknown ion case 2");
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 --> we do not extrapolate to low energies, we extrapolate
120 // to high energies (??)
121 {
122 // Do not sample between incident neutron energies:
123 //---------------------------------------------------------
124 // CosTheta:
125 G4ParticleHPVector theThVec;
126 theThVec.SetInterpolationManager(theSecondManager[it]);
127 for (i = 0; i < nCosTh[it]; i++) {
128 theThVec.SetX(i, theData[it][i].GetLabel());
129 theThVec.SetY(i, theData[it][i].GetIntegral());
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()) break;
141 }
142 // now get the prob at this energy for the right theta value
143 x = cosTh;
144 x1 = theData[it][i1 - 1].GetLabel();
145 x2 = theData[it][i1].GetLabel();
146 G4ParticleHPVector theBuff1a;
147 theBuff1a.SetInterpolationManager(theData[it][i1 - 1].GetInterpolationManager());
148 for (i = 0; i < theData[it][i1 - 1].GetVectorLength(); i++) {
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[it].GetScheme(i1), x, x1, x2, y1, y2);
153 theBuff1a.SetData(i, E, y);
154 }
155 G4ParticleHPVector theBuff2a;
156 theBuff2a.SetInterpolationManager(theData[it][i1].GetInterpolationManager());
157 for (i = 0; i < theData[it][i1].GetVectorLength(); i++) {
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[it].GetScheme(i1), x, x1, x2, y1, y2);
162 theBuff2a.SetData(i, E, y); // wrong E, right theta.
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 select theta.
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].GetIntegral());
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 neutron energy
185 x = anEnergy;
186 x1 = theEnergies[it - 1];
187 x2 = theEnergies[it];
188 G4ParticleHPVector thBuff1; // to be interpolated as run1.
189 thBuff1.SetInterpolationManager(theSecondManager[it - 1]);
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.GetScheme(it), x, x1, x2, y1, y2);
195 thBuff1.SetData(i, tmp, y);
196 }
197 G4ParticleHPVector thBuff2;
198 thBuff2.SetInterpolationManager(theSecondManager[it]);
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.GetScheme(it), x, x1, x2, y1, y2);
204 thBuff2.SetData(i, tmp, y);
205 }
206 G4ParticleHPVector theThVec;
207 theThVec.Merge(&thBuff1, &thBuff2); // takes care of interpolation
208 cosTh = theThVec.Sample();
209 /*
210 if(massCode>0.99 && massCode<1.01){theThVec.Dump();}
211 G4double random = (theThVec.GetY(theThVec.GetVectorLength()-1)
212 -theThVec.GetY(0)) *G4UniformRand();
213 G4cout<<" -- "<<theThVec.GetY(theThVec.GetVectorLength()-1)<<" "<<theThVec.GetY(0)<<" ---->
214 "<<random<<G4endl; G4int ith(0); for(i=1;i<theThVec.GetVectorLength(); i++)
215 {
216 ith = i;
217 if(random<theThVec.GetY(i)-theThVec.GetY(0)) break;
218 }
219 {
220 // calculate theta
221 G4double xx, xx1, xx2, yy1, yy2;
222 xx = random;
223 xx1 = theThVec.GetY(ith-1)-theThVec.GetY(0); // integrals
224 xx2 = theThVec.GetY(ith)-theThVec.GetY(0);
225 yy1 = theThVec.GetX(ith-1); // std::cos(theta)
226 yy2 = theThVec.GetX(ith);
227 cosTh = theInt.Interpolate(theSecondManager[it].GetScheme(ith),
228 xx, xx1,xx2,yy1,yy2);
229 }
230 */
231 G4int i1(0), i2(0);
232 // get the indixes of the vectors close to theta for low energy
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()) break;
237 }
238 // now get the prob at this energy for the right theta value
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[it - 1][i1 - 1].GetInterpolationManager());
244 for (i = 0; i < theData[it - 1][i1 - 1].GetVectorLength(); i++) {
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[it - 1].GetScheme(i1), x, x1, x2, y1, y2);
249 theBuff1a.SetData(i, E, y); // wrong E, right theta.
250 }
251 G4ParticleHPVector theBuff2a;
252 theBuff2a.SetInterpolationManager(theData[it - 1][i1].GetInterpolationManager());
253 for (i = 0; i < theData[it - 1][i1].GetVectorLength(); i++) {
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[it - 1].GetScheme(i1), x, x1, x2, y1, y2);
258 theBuff2a.SetData(i, E, y); // wrong E, right theta.
259 }
260 G4ParticleHPVector theStore1;
261 theStore1.Merge(&theBuff1a, &theBuff2a); // wrong E, right theta, complete binning
262
263 // get the indixes of the vectors close to theta for high energy
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()) break;
268 } // sonderfaelle mit i1 oder i2 head on fehlen. @@@@@
269 x1 = theData[it][i2 - 1].GetLabel();
270 x2 = theData[it][i2].GetLabel();
271 G4ParticleHPVector theBuff1b;
272 theBuff1b.SetInterpolationManager(theData[it][i2 - 1].GetInterpolationManager());
273 for (i = 0; i < theData[it][i2 - 1].GetVectorLength(); i++) {
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[it].GetScheme(i2), x, x1, x2, y1, y2);
278 theBuff1b.SetData(i, E, y); // wrong E, right theta.
279 }
280 G4ParticleHPVector theBuff2b;
281 theBuff2b.SetInterpolationManager(theData[it][i2].GetInterpolationManager());
282 // 080808 i1 -> i2
283 // for(i=0;i<theData[it][i1].GetVectorLength(); i++)
284 for (i = 0; i < theData[it][i2].GetVectorLength(); i++) {
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[it].GetScheme(i2), x, x1, x2, y1, y2);
292 theBuff2b.SetData(i, E, y); // wrong E, right theta.
293 }
294 G4ParticleHPVector theStore2;
295 theStore2.Merge(&theBuff1b, &theBuff2b); // wrong E, right theta, complete binning
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.GetInterpolationManager());
303 for (i = 0; i < theStore1.GetVectorLength(); i++) {
304 E = theStore1.GetX(i);
305 y1 = theStore1.GetY(i);
306 y2 = theStore2.GetY(E);
307 y = theInt.Interpolate(theManager.GetScheme(it), x, x1, x2, y1, y2);
308 theOne1.SetData(i, E, y); // both correct
309 }
310 G4ParticleHPVector theOne2;
311 theOne2.SetInterpolationManager(theStore2.GetInterpolationManager());
312 for (i = 0; i < theStore2.GetVectorLength(); i++) {
313 E = theStore2.GetX(i);
314 y1 = theStore1.GetY(E);
315 y2 = theStore2.GetY(i);
316 y = theInt.Interpolate(theManager.GetScheme(it), x, x1, x2, y1, y2);
317 theOne2.SetData(i, E, y); // both correct
318 }
319 G4ParticleHPVector theOne;
320 theOne.Merge(&theOne1, &theOne2); // both correct, complete binning
321
322 secEnergy = theOne.Sample();
323 currentMeanEnergy = theOne.GetMeanX();
324 }
325
326 // now do random direction in phi, and fill the result.
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::cos(phi), mtot * sinth * std::sin(phi),
335 mtot * std::cos(theta));
336 result->SetMomentum(tempVector);
337
338 return result;
339 }
340