source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPDiscreteTwoBody.cc@ 1350

Last change on this file since 1350 was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 11.7 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//080612 Bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #2,3
31//080709 Bug fix Sampling Legendre expansion by T. Koi
32//101110 Bug fix in MF=6, LAW=2 case; contribution from E. Mendoza, D. Cano-Ott (CIEMAT)
33//
34#include "G4NeutronHPDiscreteTwoBody.hh"
35#include "G4Gamma.hh"
36#include "G4Electron.hh"
37#include "G4Positron.hh"
38#include "G4Neutron.hh"
39#include "G4Proton.hh"
40#include "G4Deuteron.hh"
41#include "G4Triton.hh"
42#include "G4He3.hh"
43#include "G4Alpha.hh"
44#include "G4NeutronHPVector.hh"
45#include "G4NeutronHPLegendreStore.hh"
46
47G4ReactionProduct * G4NeutronHPDiscreteTwoBody::Sample(G4double anEnergy, G4double massCode, G4double )
48{ // Interpolation still only for the most used parts; rest to be Done @@@@@
49 G4ReactionProduct * result = new G4ReactionProduct;
50 G4int Z = static_cast<G4int>(massCode/1000);
51 G4int A = static_cast<G4int>(massCode-1000*Z);
52
53 if(massCode==0)
54 {
55 result->SetDefinition(G4Gamma::Gamma());
56 }
57 else if(A==0)
58 {
59 result->SetDefinition(G4Electron::Electron());
60 if(Z==1) result->SetDefinition(G4Positron::Positron());
61 }
62 else if(A==1)
63 {
64 result->SetDefinition(G4Neutron::Neutron());
65 if(Z==1) result->SetDefinition(G4Proton::Proton());
66 }
67 else if(A==2)
68 {
69 result->SetDefinition(G4Deuteron::Deuteron());
70 }
71 else if(A==3)
72 {
73 result->SetDefinition(G4Triton::Triton());
74 if(Z==2) result->SetDefinition(G4He3::He3());
75 }
76 else if(A==4)
77 {
78 result->SetDefinition(G4Alpha::Alpha());
79 if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "Unknown ion case 1");
80 }
81 else
82 {
83 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPDiscreteTwoBody: Unknown ion case 2");
84 }
85
86// get cosine(theta)
87 G4int i(0), it(0);
88 G4double cosTh(0);
89 for(i=0; i<nEnergy; i++)
90 {
91 it = i;
92 if(theCoeff[i].GetEnergy()>anEnergy) break;
93 }
94 if(it==0||it==nEnergy-1)
95 {
96 if(theCoeff[it].GetRepresentation()==0)
97 {
98//TK Legendre expansion
99 G4NeutronHPLegendreStore theStore(1);
100 theStore.SetCoeff(0, theCoeff);
101 theStore.SetManager(theManager);
102 //cosTh = theStore.SampleMax(anEnergy);
103 //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #3
104 cosTh = theStore.SampleDiscreteTwoBody(anEnergy);
105 }
106 else if(theCoeff[it].GetRepresentation()==12) // means LINLIN
107 {
108 G4NeutronHPVector theStore;
109 G4InterpolationManager aManager;
110 aManager.Init(LINLIN, theCoeff[it].GetNumberOfPoly()/2);
111 theStore.SetInterpolationManager(aManager);
112 for(i=0;i<theCoeff[it].GetNumberOfPoly(); i++)
113 {
114 //101110
115 //theStore.SetX(i, theCoeff[it].GetCoeff(i));
116 //theStore.SetY(i, theCoeff[it].GetCoeff(i));
117 theStore.SetX(i/2, theCoeff[it].GetCoeff(i));
118 theStore.SetY(i/2, theCoeff[it].GetCoeff(i+1));
119 i++;
120 }
121 cosTh = theStore.Sample();
122 }
123 else if(theCoeff[it].GetRepresentation()==14) //this is LOGLIN
124 {
125 G4NeutronHPVector theStore;
126 G4InterpolationManager aManager;
127 aManager.Init(LOGLIN, theCoeff[it].GetNumberOfPoly()/2);
128 theStore.SetInterpolationManager(aManager);
129 for(i=0;i<theCoeff[it].GetNumberOfPoly(); i++)
130 {
131 //101110
132 //theStore.SetX(i, theCoeff[it].GetCoeff(i));
133 //theStore.SetY(i, theCoeff[it].GetCoeff(i));
134 theStore.SetX(i/2, theCoeff[it].GetCoeff(i));
135 theStore.SetY(i/2, theCoeff[it].GetCoeff(i+1));
136 i++;
137 }
138 cosTh = theStore.Sample();
139 }
140 else
141 {
142 throw G4HadronicException(__FILE__, __LINE__, "unknown representation type in Two-body scattering");
143 }
144 }
145 else
146 {
147 if(theCoeff[it].GetRepresentation() == theCoeff[it-1].GetRepresentation())
148 {
149 if(theCoeff[it].GetRepresentation()==0)
150 {
151//TK Legendre expansion
152 G4NeutronHPLegendreStore theStore(2);
153 theStore.SetCoeff(0, &(theCoeff[it-1]));
154 theStore.SetCoeff(1, &(theCoeff[it]));
155 G4InterpolationManager aManager;
156 aManager.Init(theManager.GetScheme(it), 2);
157 theStore.SetManager(aManager);
158 //cosTh = theStore.SampleMax(anEnergy);
159//080709 TKDB
160 cosTh = theStore.SampleDiscreteTwoBody(anEnergy);
161 }
162 else if(theCoeff[it].GetRepresentation()==12) // LINLIN
163 {
164 G4NeutronHPVector theBuff1;
165 G4InterpolationManager aManager1;
166 aManager1.Init(LINLIN, theCoeff[it-1].GetNumberOfPoly()/2);
167 theBuff1.SetInterpolationManager(aManager1);
168 for(i=0;i<theCoeff[it-1].GetNumberOfPoly(); i++)
169 {
170 //101110
171 //theBuff1.SetX(i, theCoeff[it-1].GetCoeff(i));
172 //theBuff1.SetY(i, theCoeff[it-1].GetCoeff(i));
173 theBuff1.SetX(i/2, theCoeff[it-1].GetCoeff(i));
174 theBuff1.SetY(i/2, theCoeff[it-1].GetCoeff(i+1));
175 i++;
176 }
177 G4NeutronHPVector theBuff2;
178 G4InterpolationManager aManager2;
179 aManager2.Init(LINLIN, theCoeff[it].GetNumberOfPoly()/2);
180 theBuff2.SetInterpolationManager(aManager2);
181 for(i=0;i<theCoeff[it].GetNumberOfPoly(); i++)
182 {
183 //theBuff2.SetX(i, theCoeff[it].GetCoeff(i));
184 //theBuff2.SetY(i, theCoeff[it].GetCoeff(i));
185 theBuff2.SetX(i, theCoeff[it].GetCoeff(i));
186 theBuff2.SetY(i, theCoeff[it].GetCoeff(i+1));
187 i++;
188 }
189
190 G4double x1 = theCoeff[it-1].GetEnergy();
191 G4double x2 = theCoeff[it].GetEnergy();
192 G4double x = anEnergy;
193 G4double y1, y2, y, mu;
194
195 G4NeutronHPVector theStore1;
196 theStore1.SetInterpolationManager(aManager1);
197 G4NeutronHPVector theStore2;
198 theStore2.SetInterpolationManager(aManager2);
199 G4NeutronHPVector theStore;
200
201 // for fixed mu get p1, p2 and interpolate according to x
202 for(i=0; i<theBuff1.GetVectorLength(); i++)
203 {
204 mu = theBuff1.GetX(i);
205 y1 = theBuff1.GetY(i);
206 y2 = theBuff2.GetY(mu);
207 y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
208 theStore1.SetData(i, mu, y);
209 }
210 for(i=0; i<theBuff2.GetVectorLength(); i++)
211 {
212 mu = theBuff2.GetX(i);
213 y1 = theBuff2.GetY(i);
214 y2 = theBuff1.GetY(mu);
215 y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
216 theStore2.SetData(i, mu, y);
217 }
218 theStore.Merge(&theStore1, &theStore2); // merge takes care of interpolationschemes
219 cosTh = theStore.Sample();
220 }
221 else if(theCoeff[it].GetRepresentation()==14) //TK LOG_LIN
222 {
223 G4NeutronHPVector theBuff1;
224 G4InterpolationManager aManager1;
225 aManager1.Init(LOGLIN, theCoeff[it-1].GetNumberOfPoly()/2);
226 theBuff1.SetInterpolationManager(aManager1);
227 for(i=0;i<theCoeff[it-1].GetNumberOfPoly(); i++)
228 {
229 //101110
230 //theBuff1.SetX(i, theCoeff[it-1].GetCoeff(i));
231 //theBuff1.SetY(i, theCoeff[it-1].GetCoeff(i));
232 theBuff1.SetX(i/2, theCoeff[it-1].GetCoeff(i));
233 theBuff1.SetY(i/2, theCoeff[it-1].GetCoeff(i+1));
234 i++;
235 }
236
237 G4NeutronHPVector theBuff2;
238 G4InterpolationManager aManager2;
239 aManager2.Init(LOGLIN, theCoeff[it].GetNumberOfPoly()/2);
240 theBuff2.SetInterpolationManager(aManager2);
241 for(i=0;i<theCoeff[it].GetNumberOfPoly(); i++)
242 {
243 //101110
244 //theBuff2.SetX(i, theCoeff[it].GetCoeff(i));
245 //theBuff2.SetY(i, theCoeff[it].GetCoeff(i));
246 theBuff2.SetX(i/2, theCoeff[it].GetCoeff(i));
247 theBuff2.SetY(i/2, theCoeff[it].GetCoeff(i+1));
248 i++;
249 }
250
251 G4double x1 = theCoeff[it-1].GetEnergy();
252 G4double x2 = theCoeff[it].GetEnergy();
253 G4double x = anEnergy;
254 G4double y1, y2, y, mu;
255
256 G4NeutronHPVector theStore1;
257 theStore1.SetInterpolationManager(aManager1);
258 G4NeutronHPVector theStore2;
259 theStore2.SetInterpolationManager(aManager2);
260 G4NeutronHPVector theStore;
261
262 // for fixed mu get p1, p2 and interpolate according to x
263 for(i=0; i<theBuff1.GetVectorLength(); i++)
264 {
265 mu = theBuff1.GetX(i);
266 y1 = theBuff1.GetY(i);
267 y2 = theBuff2.GetY(mu);
268 y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
269 theStore1.SetData(i, mu, y);
270 }
271 for(i=0; i<theBuff2.GetVectorLength(); i++)
272 {
273 mu = theBuff2.GetX(i);
274 y1 = theBuff2.GetY(i);
275 y2 = theBuff1.GetY(mu);
276 y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
277 theStore2.SetData(i, mu, y);
278 }
279 theStore.Merge(&theStore1, &theStore2);
280 cosTh = theStore.Sample();
281 }
282 else
283 {
284 throw G4HadronicException(__FILE__, __LINE__, "Two neighbouring distributions with different interpolation");
285 }
286 }
287 else
288 {
289 throw G4HadronicException(__FILE__, __LINE__, "unknown representation type in Two-body scattering, case 2");
290 }
291 }
292
293// now get the energy from kinematics and Q-value.
294
295 //G4double restEnergy = anEnergy+GetQValue();
296
297// assumed to be in CMS @@@@@@@@@@@@@@@@@
298
299 //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #2
300 //G4double residualMass = GetTarget()->GetMass() + GetNeutron()->GetMass()
301 // - result->GetMass() - GetQValue();
302 //G4double kinE = restEnergy/(1+result->GetMass()/residualMass); // non relativistic @@
303 G4double A1 = GetTarget()->GetMass()/GetNeutron()->GetMass();
304 G4double A1prim = result->GetMass()/GetNeutron()->GetMass();
305 G4double E1 = (A1+1)*(A1+1)/A1/A1*anEnergy;
306 G4double kinE = (A1+1-A1prim)/(A1+1)/(A1+1)*(A1*E1+(1+A1)*GetQValue());
307
308 result->SetKineticEnergy(kinE); // non relativistic @@
309 G4double phi = twopi*G4UniformRand();
310 G4double theta = std::acos(cosTh);
311 G4double sinth = std::sin(theta);
312 G4double mtot = result->GetTotalMomentum();
313 G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
314 result->SetMomentum(tempVector);
315
316// some garbage collection
317
318// return the result
319 return result;
320}
Note: See TracBrowser for help on using the repository browser.