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