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

Last change on this file since 1347 was 1347, checked in by garnier, 13 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.