source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPLabAngularEnergy.cc@ 893

Last change on this file since 893 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 12.6 KB
Line 
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#include "G4NeutronHPLabAngularEnergy.hh"
31#include "G4Gamma.hh"
32#include "G4Electron.hh"
33#include "G4Positron.hh"
34#include "G4Neutron.hh"
35#include "G4Proton.hh"
36#include "G4Deuteron.hh"
37#include "G4Triton.hh"
38#include "G4He3.hh"
39#include "G4Alpha.hh"
40#include "Randomize.hh"
41
42void G4NeutronHPLabAngularEnergy::Init(std::ifstream & aDataFile)
43{
44 aDataFile >> nEnergies;
45 theManager.Init(aDataFile);
46 theEnergies = new G4double[nEnergies];
47 nCosTh = new G4int[nEnergies];
48 theData = new G4NeutronHPVector * [nEnergies];
49 theSecondManager = new G4InterpolationManager [nEnergies];
50 for(G4int i=0; i<nEnergies; i++)
51 {
52 aDataFile >> theEnergies[i];
53 theEnergies[i]*=eV;
54 aDataFile >> nCosTh[i];
55 theSecondManager[i].Init(aDataFile);
56 theData[i] = new G4NeutronHPVector[nCosTh[i]];
57 G4double label;
58 for(G4int ii=0; ii<nCosTh[i]; ii++)
59 {
60 aDataFile >> label;
61 theData[i][ii].SetLabel(label);
62 theData[i][ii].Init(aDataFile, eV);
63 }
64 }
65}
66
67G4ReactionProduct * G4NeutronHPLabAngularEnergy::Sample(G4double anEnergy, G4double massCode, G4double )
68{
69 G4ReactionProduct * result = new G4ReactionProduct;
70 G4int Z = static_cast<G4int>(massCode/1000);
71 G4int A = static_cast<G4int>(massCode-1000*Z);
72
73 if(massCode==0)
74 {
75 result->SetDefinition(G4Gamma::Gamma());
76 }
77 else if(A==0)
78 {
79 result->SetDefinition(G4Electron::Electron());
80 if(Z==1) result->SetDefinition(G4Positron::Positron());
81 }
82 else if(A==1)
83 {
84 result->SetDefinition(G4Neutron::Neutron());
85 if(Z==1) result->SetDefinition(G4Proton::Proton());
86 }
87 else if(A==2)
88 {
89 result->SetDefinition(G4Deuteron::Deuteron());
90 }
91 else if(A==3)
92 {
93 result->SetDefinition(G4Triton::Triton());
94 if(Z==2) result->SetDefinition(G4He3::He3());
95 }
96 else if(A==4)
97 {
98 result->SetDefinition(G4Alpha::Alpha());
99 if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "Unknown ion case 1");
100 }
101 else
102 {
103 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPLabAngularEnergy: Unknown ion case 2");
104 }
105
106 // get theta, E
107 G4double cosTh, secEnergy;
108 G4int i, it(0);
109 // find the energy bin
110 for(i=0; i<nEnergies; i++)
111 {
112 it = i;
113 if(anEnergy<theEnergies[i]) break;
114 }
115 if(it==0 || it == nEnergies-1) // it marks the energy bin
116 {
117 // integrate the prob for each costh, and select theta.
118 G4double * running = new G4double [nCosTh[it]];
119 running[0]=0;
120 for(i=0;i<nCosTh[it]; i++)
121 {
122 if(i!=0) running[i] = running[i-1];
123 running[i]+=theData[it][i].GetIntegral(); // Does interpolated integral.
124 }
125 G4double random = running[nCosTh[it]-1]*G4UniformRand();
126 G4int ith(0);
127 for(i=0;i<nCosTh[it]; i++)
128 {
129 ith = i;
130 if(random<running[i]) break;
131 }
132 if(ith==0 || ith==nCosTh[it]-1)
133 {
134 cosTh = theData[it][ith].GetLabel();
135 secEnergy = theData[it][ith].Sample();
136 currentMeanEnergy = theData[it][ith].GetMeanX();
137 }
138 else
139 {
140 G4double x1 = theData[it][ith-1].GetIntegral();
141 G4double x2 = theData[it][ith].GetIntegral();
142 G4double x = random;
143 G4double y1 = theData[it][ith-1].GetLabel();
144 G4double y2 = theData[it][ith].GetLabel();
145 cosTh = theInt.Interpolate(theSecondManager[it].GetInverseScheme(ith),
146 x, x1, x2, y1, y2);
147 G4NeutronHPVector theBuff1;
148 theBuff1.SetInterpolationManager(theData[it][ith-1].GetInterpolationManager());
149 G4NeutronHPVector theBuff2;
150 theBuff2.SetInterpolationManager(theData[it][ith].GetInterpolationManager());
151 x1=y1;
152 x2=y2;
153 G4double y, mu;
154 for(i=0;i<theData[it][ith-1].GetVectorLength(); i++)
155 {
156 mu = theData[it][ith-1].GetX(i);
157 y1 = theData[it][ith-1].GetY(i);
158 y2 = theData[it][ith].GetY(mu);
159 y = theInt.Interpolate(theSecondManager[it].GetScheme(ith),
160 cosTh, x1,x2,y1,y2);
161 theBuff1.SetData(i, mu, y);
162 }
163 for(i=0;i<theData[it][ith].GetVectorLength(); i++)
164 {
165 mu = theData[it][ith].GetX(i);
166 y1 = theData[it][ith-1].GetY(mu);
167 y2 = theData[it][ith].GetY(i);
168 y = theInt.Interpolate(theSecondManager[it].GetScheme(ith),
169 cosTh, x1,x2,y1,y2);
170 theBuff2.SetData(i, mu, y);
171 }
172 G4NeutronHPVector theStore;
173 theStore.Merge(&theBuff1, &theBuff2);
174 secEnergy = theStore.Sample();
175 currentMeanEnergy = theStore.GetMeanX();
176 }
177 delete [] running;
178 }
179 else // this is the small big else.
180 {
181 G4double x, x1, x2, y1, y2, y, tmp, E;
182 // integrate the prob for each costh, and select theta.
183 G4NeutronHPVector run1;
184 run1.SetY(0, 0.);
185 for(i=0;i<nCosTh[it-1]; i++)
186 {
187 if(i!=0) run1.SetY(i, run1.GetY(i-1));
188 run1.SetX(i, theData[it-1][i].GetLabel());
189 run1.SetY(i, run1.GetY(i)+theData[it-1][i].GetIntegral());
190 }
191 G4NeutronHPVector run2;
192 run2.SetY(0, 0.);
193 for(i=0;i<nCosTh[it]; i++)
194 {
195 if(i!=0) run2.SetY(i, run2.GetY(i-1));
196 run2.SetX(i, theData[it][i].GetLabel());
197 run2.SetY(i, run2.GetY(i)+theData[it][i].GetIntegral());
198 }
199 // get the distributions for the correct neutron energy
200 x = anEnergy;
201 x1 = theEnergies[it-1];
202 x2 = theEnergies[it];
203 G4NeutronHPVector thBuff1; // to be interpolated as run1.
204 thBuff1.SetInterpolationManager(theSecondManager[it-1]);
205 for(i=0; i<run1.GetVectorLength(); i++)
206 {
207 tmp = run1.GetX(i); //theta
208 y1 = run1.GetY(i); // integral
209 y2 = run2.GetY(tmp);
210 y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
211 thBuff1.SetData(i, tmp, y);
212 }
213 G4NeutronHPVector thBuff2;
214 thBuff2.SetInterpolationManager(theSecondManager[it]);
215 for(i=0; i<run2.GetVectorLength(); i++)
216 {
217 tmp = run2.GetX(i); //theta
218 y1 = run1.GetY(tmp); // integral
219 y2 = run2.GetY(i);
220 y = theInt.Lin(x, x1,x2,y1,y2);
221 thBuff2.SetData(i, tmp, y);
222 }
223 G4NeutronHPVector theThVec;
224 theThVec.Merge(&thBuff1 ,&thBuff2); // takes care of interpolation
225 G4double random = (theThVec.GetY(theThVec.GetVectorLength()-1)
226 -theThVec.GetY(0)) *G4UniformRand();
227 G4int ith(0);
228 for(i=1;i<theThVec.GetVectorLength(); i++)
229 {
230 ith = i;
231 if(random<theThVec.GetY(i)-theThVec.GetY(0)) break;
232 }
233 {
234 // calculate theta
235 G4double x, x1, x2, y1, y2;
236 x = random;
237 x1 = theThVec.GetY(ith-1)-theThVec.GetY(0); // integrals
238 x2 = theThVec.GetY(ith)-theThVec.GetY(0);
239 y1 = theThVec.GetX(ith-1); // std::cos(theta)
240 y2 = theThVec.GetX(ith);
241 cosTh = theInt.Interpolate(theSecondManager[it].GetScheme(ith),
242 x, x1,x2,y1,y2);
243 }
244 G4int i1(0), i2(0);
245 // get the indixes of the vectors close to theta for low energy
246 // first it-1 !!!! i.e. low in energy
247 for(i=0; i<nCosTh[it-1]; i++)
248 {
249 i1 = i;
250 if(cosTh<theData[it-1][i].GetLabel()) break;
251 }
252 // now get the prob at this energy for the right theta value
253 x = cosTh;
254 x1 = theData[it-1][i1-1].GetLabel();
255 x2 = theData[it-1][i1].GetLabel();
256 G4NeutronHPVector theBuff1a;
257 theBuff1a.SetInterpolationManager(theData[it-1][i1-1].GetInterpolationManager());
258 for(i=0;i<theData[it-1][i1-1].GetVectorLength(); i++)
259 {
260 E = theData[it-1][i1-1].GetX(i);
261 y1 = theData[it-1][i1-1].GetY(i);
262 y2 = theData[it-1][i1].GetY(E);
263 y = theInt.Lin(x, x1,x2,y1,y2);
264 theBuff1a.SetData(i, E, y); // wrong E, right theta.
265 }
266 G4NeutronHPVector theBuff2a;
267 theBuff2a.SetInterpolationManager(theData[it-1][i1].GetInterpolationManager());
268 for(i=0;i<theData[it-1][i1].GetVectorLength(); i++)
269 {
270 E = theData[it-1][i1].GetX(i);
271 y1 = theData[it-1][i1-1].GetY(E);
272 y2 = theData[it-1][i1].GetY(i);
273 y = theInt.Lin(x, x1,x2,y1,y2);
274 theBuff2a.SetData(i, E, y); // wrong E, right theta.
275 }
276 G4NeutronHPVector theStore1;
277 theStore1.Merge(&theBuff1a, &theBuff2a); // wrong E, right theta, complete binning
278
279 // get the indixes of the vectors close to theta for high energy
280 // then it !!!! i.e. high in energy
281 for(i=0; i<nCosTh[it]; i++)
282 {
283 i2 = i;
284 if(cosTh<theData[it][i2].GetLabel()) break;
285 } // sonderfaelle mit i1 oder i2 head on fehlen. @@@@@
286 x1 = theData[it][i2-1].GetLabel();
287 x2 = theData[it][i2].GetLabel();
288 G4NeutronHPVector theBuff1b;
289 theBuff1b.SetInterpolationManager(theData[it][i2-1].GetInterpolationManager());
290 for(i=0;i<theData[it][i2-1].GetVectorLength(); i++)
291 {
292 E = theData[it][i2-1].GetX(i);
293 y1 = theData[it][i2-1].GetY(i);
294 y2 = theData[it][i2].GetY(E);
295 y = theInt.Lin(x, x1,x2,y1,y2);
296 theBuff1b.SetData(i, E, y); // wrong E, right theta.
297 }
298 G4NeutronHPVector theBuff2b;
299 theBuff2b.SetInterpolationManager(theData[it][i2].GetInterpolationManager());
300 for(i=0;i<theData[it][i1].GetVectorLength(); i++)
301 {
302 E = theData[it][i1].GetX(i);
303 y1 = theData[it][i1-1].GetY(E);
304 y2 = theData[it][i1].GetY(i);
305 y = theInt.Lin(x, x1,x2,y1,y2);
306 theBuff2b.SetData(i, E, y); // wrong E, right theta.
307 }
308 G4NeutronHPVector theStore2;
309 theStore2.Merge(&theBuff1b, &theBuff2b); // wrong E, right theta, complete binning
310 // now get to the right energy.
311
312 x = anEnergy;
313 x1 = theEnergies[it-1];
314 x2 = theEnergies[it];
315 G4NeutronHPVector theOne1;
316 theOne1.SetInterpolationManager(theStore1.GetInterpolationManager());
317 for(i=0; i<theStore1.GetVectorLength(); i++)
318 {
319 E = theStore1.GetX(i);
320 y1 = theStore1.GetY(i);
321 y2 = theStore2.GetY(E);
322 y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
323 theOne1.SetData(i, E, y); // both correct
324 }
325 G4NeutronHPVector theOne2;
326 theOne2.SetInterpolationManager(theStore2.GetInterpolationManager());
327 for(i=0; i<theStore2.GetVectorLength(); i++)
328 {
329 E = theStore2.GetX(i);
330 y1 = theStore1.GetY(E);
331 y2 = theStore2.GetY(i);
332 y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
333 theOne2.SetData(i, E, y); // both correct
334 }
335 G4NeutronHPVector theOne;
336 theOne.Merge(&theOne1, &theOne2); // both correct, complete binning
337
338 secEnergy = theOne.Sample();
339 currentMeanEnergy = theOne.GetMeanX();
340 }
341
342// now do random direction in phi, and fill the result.
343
344 result->SetKineticEnergy(secEnergy);
345
346 G4double phi = twopi*G4UniformRand();
347 G4double theta = std::acos(cosTh);
348 G4double sinth = std::sin(theta);
349 G4double mtot = result->GetTotalMomentum();
350 G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
351 result->SetMomentum(tempVector);
352
353 return result;
354}
Note: See TracBrowser for help on using the repository browser.