Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
G4NeutronHPDiscreteTwoBody Class Reference

#include <G4NeutronHPDiscreteTwoBody.hh>

Inheritance diagram for G4NeutronHPDiscreteTwoBody:
G4VNeutronHPEnergyAngular

Public Member Functions

 G4NeutronHPDiscreteTwoBody ()
 
 ~G4NeutronHPDiscreteTwoBody ()
 
void Init (std::istream &aDataFile)
 
G4ReactionProduct * Sample (G4double anEnergy, G4double massCode, G4double mass)
 
G4double MeanEnergyOfThisInteraction ()
 
- Public Member Functions inherited from G4VNeutronHPEnergyAngular
 G4VNeutronHPEnergyAngular ()
 
virtual ~G4VNeutronHPEnergyAngular ()
 
void SetNeutron (G4ReactionProduct *aNeutron)
 
void SetTarget (G4ReactionProduct *aTarget)
 
G4ReactionProduct * GetTarget ()
 
G4ReactionProduct * GetNeutron ()
 
G4ReactionProduct * GetCMS ()
 
void SetQValue (G4double aValue)
 
virtual void ClearHistories ()
 

Additional Inherited Members

- Protected Member Functions inherited from G4VNeutronHPEnergyAngular
G4double GetQValue ()
 

Detailed Description

Definition at line 41 of file G4NeutronHPDiscreteTwoBody.hh.

Constructor & Destructor Documentation

G4NeutronHPDiscreteTwoBody::G4NeutronHPDiscreteTwoBody ( )
inline

Definition at line 45 of file G4NeutronHPDiscreteTwoBody.hh.

46  {
47  theCoeff = 0;
48  }
G4NeutronHPDiscreteTwoBody::~G4NeutronHPDiscreteTwoBody ( )
inline

Definition at line 49 of file G4NeutronHPDiscreteTwoBody.hh.

50  {
51  if(theCoeff!=0) delete [] theCoeff;
52  }

Member Function Documentation

void G4NeutronHPDiscreteTwoBody::Init ( std::istream &  aDataFile)
inlinevirtual

Implements G4VNeutronHPEnergyAngular.

Definition at line 54 of file G4NeutronHPDiscreteTwoBody.hh.

References energy(), G4NeutronHPLegendreTable::Init(), G4InterpolationManager::Init(), G4NeutronHPLegendreTable::SetCoeff(), and G4NeutronHPLegendreTable::SetRepresentation().

55  {
56  aDataFile >> nEnergy;
57  theManager.Init(aDataFile);
58  theCoeff = new G4NeutronHPLegendreTable[nEnergy];
59  for(G4int i=0; i<nEnergy; i++)
60  {
61  G4double energy;
62  G4int aRep, nCoeff;
63  aDataFile >> energy >> aRep >> nCoeff;
64  energy*=CLHEP::eV;
65  G4int nPoints=nCoeff;
66  if(aRep>0) nPoints*=2;
67  //theCoeff[i].Init(energy, nPoints);
68 
69  theCoeff[i].Init(energy, nPoints-1);
70  theCoeff[i].SetRepresentation(aRep);
71  for(G4int ii=0; ii<nPoints; ii++)
72  {
73  G4double y;
74  aDataFile >> y;
75  theCoeff[i].SetCoeff(ii, y);
76  }
77  }
78  }
G4InterpolationManager::Init
void Init(G4int aScheme, G4int aRange)
Definition: G4InterpolationManager.hh:84
G4NeutronHPLegendreTable::Init
void Init(std::istream &aDataFile)
Definition: G4NeutronHPLegendreTable.hh:66
G4int
int G4int
Definition: G4Types.hh:78
energy
double precision function energy(A, Z)
Definition: dpm25nuc6.f:4106
G4NeutronHPLegendreTable::SetCoeff
void SetCoeff(G4int l, G4double coeff)
Definition: G4NeutronHPLegendreTable.hh:91
G4double
double G4double
Definition: G4Types.hh:76
G4double G4NeutronHPDiscreteTwoBody::MeanEnergyOfThisInteraction ( )
inlinevirtual

Implements G4VNeutronHPEnergyAngular.

Definition at line 81 of file G4NeutronHPDiscreteTwoBody.hh.

81 { return -1; }
G4ReactionProduct * G4NeutronHPDiscreteTwoBody::Sample ( G4double  anEnergy,
G4double  massCode,
G4double  mass 
)
virtual

Implements G4VNeutronHPEnergyAngular.

Definition at line 48 of file G4NeutronHPDiscreteTwoBody.cc.

References G4Alpha::Alpha(), G4Deuteron::Deuteron(), G4Electron::Electron(), G4UniformRand, G4Gamma::Gamma(), G4NeutronHPLegendreTable::GetEnergy(), G4ReactionProduct::GetMass(), G4VNeutronHPEnergyAngular::GetNeutron(), G4NeutronHPLegendreTable::GetNumberOfPoly(), G4VNeutronHPEnergyAngular::GetQValue(), G4InterpolationManager::GetScheme(), G4VNeutronHPEnergyAngular::GetTarget(), G4ReactionProduct::GetTotalMomentum(), G4NeutronHPVector::GetVectorLength(), G4NeutronHPVector::GetX(), G4NeutronHPVector::GetY(), G4He3::He3(), G4InterpolationManager::Init(), G4NeutronHPInterpolator::Interpolate(), LINLIN, LOGLIN, G4NeutronHPVector::Merge(), G4Neutron::Neutron(), G4Positron::Positron(), G4Proton::Proton(), G4NeutronHPVector::Sample(), G4NeutronHPLegendreStore::SampleDiscreteTwoBody(), G4NeutronHPLegendreStore::SetCoeff(), G4NeutronHPVector::SetData(), G4ReactionProduct::SetDefinition(), G4NeutronHPVector::SetInterpolationManager(), G4ReactionProduct::SetKineticEnergy(), G4NeutronHPLegendreStore::SetManager(), G4ReactionProduct::SetMomentum(), G4NeutronHPVector::SetX(), G4NeutronHPVector::SetY(), G4Triton::Triton(), python.hepunit::twopi, and test::x.

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 }
G4NeutronHPVector::GetY
G4double GetY(G4double x)
Definition: G4NeutronHPVector.hh:192
G4NeutronHPVector::GetVectorLength
G4int GetVectorLength() const
Definition: G4NeutronHPVector.hh:193
G4NeutronHPVector::SetInterpolationManager
void SetInterpolationManager(const G4InterpolationManager &aManager)
Definition: G4NeutronHPVector.hh:477
G4ReactionProduct::GetTotalMomentum
G4double GetTotalMomentum() const
Definition: G4ReactionProduct.hh:126
G4InterpolationManager::Init
void Init(G4int aScheme, G4int aRange)
Definition: G4InterpolationManager.hh:84
G4ReactionProduct::SetKineticEnergy
void SetKineticEnergy(const G4double en)
Definition: G4ReactionProduct.hh:132
G4ReactionProduct::SetMomentum
void SetMomentum(const G4double x, const G4double y, const G4double z)
Definition: G4ReactionProduct.cc:165
LINLIN
Definition: G4InterpolationScheme.hh:39
G4NeutronHPVector::Merge
void Merge(G4NeutronHPVector *active, G4NeutronHPVector *passive)
Definition: G4NeutronHPVector.hh:267
G4NeutronHPVector::GetX
G4double GetX(G4int i) const
Definition: G4NeutronHPVector.hh:121
G4NeutronHPVector::SetData
void SetData(G4int i, G4double x, G4double y)
Definition: G4NeutronHPVector.hh:90
G4int
int G4int
Definition: G4Types.hh:78
G4NeutronHPInterpolator::Interpolate
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
Definition: G4NeutronHPInterpolator.hh:85
G4NeutronHPVector::SetY
void SetY(G4int i, G4double x)
Definition: G4NeutronHPVector.hh:107
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4InterpolationManager::GetScheme
G4InterpolationScheme GetScheme(G4int index) const
Definition: G4InterpolationManager.hh:131
G4Triton::Triton
static G4Triton * Triton()
Definition: G4Triton.cc:95
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4NeutronHPVector::SetX
void SetX(G4int i, G4double e)
Definition: G4NeutronHPVector.hh:97
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4Deuteron::Deuteron
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4Positron::Positron
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4ReactionProduct::SetDefinition
void SetDefinition(G4ParticleDefinition *aParticleDefinition)
Definition: G4ReactionProduct.cc:154
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4double
double G4double
Definition: G4Types.hh:76
LOGLIN
Definition: G4InterpolationScheme.hh:39
G4He3::He3
static G4He3 * He3()
Definition: G4He3.cc:94
G4ReactionProduct::GetMass
G4double GetMass() const
Definition: G4ReactionProduct.hh:150
The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:56:47 for Geant4.10 by   doxygen 1.8.7