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

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

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