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

Last change on this file since 1036 was 962, checked in by garnier, 17 years ago

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