Context Navigation

source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPContAngularPar.cc @ 847

Last change on this file since 847 was 819, checked in by garnier, 16 years ago
import all except CVS
File size: 15.9 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	// 09-May-06 fix in Sample by T. Koi
31	// 080318 Fix Compilation warnings - gcc-4.3.0 by T. Koi
32	// (This fix has a real effect to the code.)
33	// 080409 Fix div0 error with G4FPE by T. Koi
34	//
35
36	#include "G4NeutronHPContAngularPar.hh"
37	#include "G4NeutronHPLegendreStore.hh"
38	#include "G4Gamma.hh"
39	#include "G4Electron.hh"
40	#include "G4Positron.hh"
41	#include "G4Neutron.hh"
42	#include "G4Proton.hh"
43	#include "G4Deuteron.hh"
44	#include "G4Triton.hh"
45	#include "G4He3.hh"
46	#include "G4Alpha.hh"
47	#include "G4NeutronHPVector.hh"
48	#include "G4NucleiPropertiesTable.hh"
49	#include "G4NeutronHPKallbachMannSyst.hh"
50	#include "G4ParticleTable.hh"
51
52	void G4NeutronHPContAngularPar::Init(std::ifstream & aDataFile)
53	{
54	aDataFile >> theEnergy >> nEnergies >> nDiscreteEnergies >> nAngularParameters;
55	theEnergy *= eV;
56	theAngular = new G4NeutronHPList [nEnergies];
57	for(G4int i=0; i<nEnergies; i++)
58	{
59	G4double sEnergy;
60	aDataFile >> sEnergy;
61	sEnergy*=eV;
62	theAngular[i].SetLabel(sEnergy);
63	theAngular[i].Init(aDataFile, nAngularParameters, 1.);
64	}
65	}
66
67	G4ReactionProduct *
68	G4NeutronHPContAngularPar::Sample(G4double anEnergy, G4double massCode, G4double /targetMass/,
69	G4int angularRep, G4int /interpolE/)
70	{
71	G4ReactionProduct * result = new G4ReactionProduct;
72	G4int Z = static_cast<G4int>(massCode/1000);
73	G4int A = static_cast<G4int>(massCode-1000*Z);
74	if(massCode==0)
75	{
76	result->SetDefinition(G4Gamma::Gamma());
77	}
78	else if(A==0)
79	{
80	result->SetDefinition(G4Electron::Electron());
81	if(Z==1) result->SetDefinition(G4Positron::Positron());
82	}
83	else if(A==1)
84	{
85	result->SetDefinition(G4Neutron::Neutron());
86	if(Z==1) result->SetDefinition(G4Proton::Proton());
87	}
88	else if(A==2)
89	{
90	result->SetDefinition(G4Deuteron::Deuteron());
91	}
92	else if(A==3)
93	{
94	result->SetDefinition(G4Triton::Triton());
95	if(Z==2) result->SetDefinition(G4He3::He3());
96	}
97	else if(A==4)
98	{
99	result->SetDefinition(G4Alpha::Alpha());
100	if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPContAngularPar: Unknown ion case 1");
101	}
102	else
103	{
104	result->SetDefinition(G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z));
105	}
106	G4int i(0);
107	G4int it(0);
108	G4double fsEnergy(0);
109	G4double cosTh(0);
110	if(angularRep==1)
111	{
112	G4double random = G4UniformRand();
113	G4double * running = new G4double[nEnergies];
114	running[0]=0;
115	G4double weighted = 0;
116	for(i=1; i<nEnergies; i++)
117	{
118	if(i!=0)
119	{
120	running[i]=running[i-1];
121	}
122	running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
123	theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
124	theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
125	weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
126	theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
127	theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
128	}
129	// cash the mean energy in this distribution
130	//080409 TKDB
131	if ( nEnergies == 1 )
132	currentMeanEnergy = 0.0;
133	else
134	currentMeanEnergy = weighted/running[nEnergies-1];
135
136	//080409 TKDB
137	if ( nEnergies == 1 ) it = 0;
138	//for(i=1; i<nEnergies; i++)
139	for(i=1; i<nEnergies; i++)
140	{
141	it = i;
142	if(random<running[i]/running[nEnergies-1]) break;
143	}
144	if(it<nDiscreteEnergies\|\|it==0)
145	{
146	if(it == 0)
147	{
148	fsEnergy = theAngular[it].GetLabel();
149	G4NeutronHPLegendreStore theStore(1);
150	theStore.Init(0,fsEnergy,nAngularParameters);
151	for(i=0;i<nAngularParameters;i++)
152	{
153	theStore.SetCoeff(0,i,theAngular[it].GetValue(i));
154	}
155	// use it to sample.
156	cosTh = theStore.SampleMax(fsEnergy);
157	}
158	else
159	{
160	G4double e1, e2;
161	e1 = theAngular[it-1].GetLabel();
162	e2 = theAngular[it].GetLabel();
163	fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
164	random,
165	running[it-1]/running[nEnergies-1],
166	running[it]/running[nEnergies-1],
167	e1, e2);
168	// fill a Legendrestore
169	G4NeutronHPLegendreStore theStore(2);
170	theStore.Init(0,e1,nAngularParameters);
171	theStore.Init(1,e2,nAngularParameters);
172	for(i=0;i<nAngularParameters;i++)
173	{
174	theStore.SetCoeff(0,i,theAngular[it-1].GetValue(i));
175	theStore.SetCoeff(1,i,theAngular[it].GetValue(i));
176	}
177	// use it to sample.
178	theStore.SetManager(theManager);
179	cosTh = theStore.SampleMax(fsEnergy);
180	}
181	}
182	else // continuum contribution
183	{
184	G4double x1 = running[it-1]/running[nEnergies-1];
185	G4double x2 = running[it]/running[nEnergies-1];
186	G4double y1 = theAngular[it-1].GetLabel();
187	G4double y2 = theAngular[it].GetLabel();
188	fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
189	random,x1,x2,y1,y2);
190	G4NeutronHPLegendreStore theStore(2);
191	theStore.Init(0,y1,nAngularParameters);
192	theStore.Init(1,y2,nAngularParameters);
193	theStore.SetManager(theManager);
194	for(i=0;i<nAngularParameters;i++)
195	{
196	theStore.SetCoeff(0,i,theAngular[it-1].GetValue(i));
197	theStore.SetCoeff(1,i,theAngular[it].GetValue(i));
198	}
199	// use it to sample.
200	cosTh = theStore.SampleMax(fsEnergy);
201	}
202	delete [] running;
203	}
204	else if(angularRep==2)
205	{
206	// first get the energy (already the right for this incoming energy)
207	G4int i;
208	G4double * running = new G4double[nEnergies];
209	running[0]=0;
210	G4double weighted = 0;
211	for(i=1; i<nEnergies; i++)
212	{
213	if(i!=0) running[i]=running[i-1];
214	running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
215	theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
216	theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
217	weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
218	theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
219	theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
220	}
221	// cash the mean energy in this distribution
222	//080409 TKDB
223	//currentMeanEnergy = weighted/running[nEnergies-1];
224	if ( nEnergies == 1 )
225	currentMeanEnergy = 0.0;
226	else
227	currentMeanEnergy = weighted/running[nEnergies-1];
228
229	G4int it(0);
230	G4double randkal = G4UniformRand();
231	//080409 TKDB
232	//for(i=0; i<nEnergies; i++)
233	for(i=1; i<nEnergies; i++)
234	{
235	it = i;
236	if(randkal<running[i]/running[nEnergies-1]) break;
237	}
238
239	// interpolate the secondary energy.
240	G4double x, x1,x2,y1,y2;
241	if(it==0) it=1;
242	x = randkal*running[nEnergies-1];
243	x1 = running[it-1];
244	x2 = running[it];
245	G4double compoundFraction;
246	// interpolate energy
247	y1 = theAngular[it-1].GetLabel();
248	y2 = theAngular[it].GetLabel();
249	fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it-1),
250	x, x1,x2,y1,y2);
251	// for theta interpolate the compoundFractions
252	G4double cLow = theAngular[it-1].GetValue(1);
253	G4double cHigh = theAngular[it].GetValue(1);
254	compoundFraction = theInt.Interpolate(theManager.GetScheme(it),
255	fsEnergy, y1, y2, cLow,cHigh);
256	delete [] running;
257
258	// get cosTh
259	G4double incidentEnergy = anEnergy;
260	G4double incidentMass = G4Neutron::Neutron()->GetPDGMass();
261	G4double productEnergy = fsEnergy;
262	G4double productMass = result->GetMass();
263	G4int targetZ = G4int(theTargetCode/1000);
264	G4int targetA = G4int(theTargetCode-1000*targetZ);
265	// To correspond to natural composition (-nat-) data files.
266	if ( targetA == 0 )
267	targetA = int ( theTarget->GetMass()/amu_c2 + 0.5 );
268	G4double targetMass = theTarget->GetMass();
269	G4int residualA = targetA+1-A;
270	G4int residualZ = targetZ-Z;
271	G4double residualMass = residualZ*G4Proton::Proton()->GetPDGMass();
272	residualMass +=(residualA-residualZ)*G4Neutron::Neutron()->GetPDGMass();
273	residualMass -= G4NucleiPropertiesTable::GetBindingEnergy(residualZ, residualA);
274	G4NeutronHPKallbachMannSyst theKallbach(compoundFraction,
275	incidentEnergy, incidentMass,
276	productEnergy, productMass,
277	residualMass, residualA, residualZ,
278	targetMass, targetA, targetZ);
279	cosTh = theKallbach.Sample(anEnergy);
280	}
281	else if(angularRep>10&&angularRep<16)
282	{
283	G4double random = G4UniformRand();
284	G4double * running = new G4double[nEnergies];
285	running[0]=0;
286	G4double weighted = 0;
287	for(i=1; i<nEnergies; i++)
288	{
289	if(i!=0) running[i]=running[i-1];
290	running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
291	theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
292	theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
293	weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
294	theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
295	theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
296	}
297	// cash the mean energy in this distribution
298	//currentMeanEnergy = weighted/running[nEnergies-1];
299	if ( nEnergies == 1 )
300	currentMeanEnergy = 0.0;
301	else
302	currentMeanEnergy = weighted/running[nEnergies-1];
303
304	//080409 TKDB
305	if ( nEnergies == 1 ) it = 0;
306	//for(i=0; i<nEnergies; i++)
307	for(i=1; i<nEnergies; i++)
308	{
309	it = i;
310	if(random<running[i]/running[nEnergies-1]) break;
311	}
312	if(it<nDiscreteEnergies\|\|it==0)
313	{
314	if(it==0)
315	{
316	fsEnergy = theAngular[0].GetLabel();
317	G4NeutronHPVector theStore;
318	G4int aCounter = 0;
319	for(G4int i=1; i<nAngularParameters; i+=2)
320	{
321	theStore.SetX(aCounter, theAngular[0].GetValue(i));
322	theStore.SetY(aCounter, theAngular[0].GetValue(i+1));
323	aCounter++;
324	}
325	G4InterpolationManager aMan;
326	aMan.Init(angularRep-10, nAngularParameters-1);
327	theStore.SetInterpolationManager(aMan);
328	cosTh = theStore.Sample();
329	}
330	else
331	{
332	fsEnergy = theAngular[it].GetLabel();
333	G4NeutronHPVector theStore;
334	G4InterpolationManager aMan;
335	aMan.Init(angularRep-10, nAngularParameters-1);
336	theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
337	G4InterpolationScheme currentScheme = theManager.GetInverseScheme(it);
338	G4int aCounter = 0;
339	for(G4int i=1; i<nAngularParameters; i+=2)
340	{
341	theStore.SetX(aCounter, theAngular[it].GetValue(i));
342	theStore.SetY(aCounter, theInt.Interpolate(currentScheme,
343	random,
344	running[it-1]/running[nEnergies-1],
345	running[it]/running[nEnergies-1],
346	theAngular[it-1].GetValue(i+1),
347	theAngular[it].GetValue(i+1)));
348	aCounter++;
349	}
350	cosTh = theStore.Sample();
351	}
352	}
353	else
354	{
355	G4double x1 = running[it-1]/running[nEnergies-1];
356	G4double x2 = running[it]/running[nEnergies-1];
357	G4double y1 = theAngular[it-1].GetLabel();
358	G4double y2 = theAngular[it].GetLabel();
359	fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
360	random,x1,x2,y1,y2);
361	G4NeutronHPVector theBuff1;
362	G4NeutronHPVector theBuff2;
363	G4InterpolationManager aMan;
364	aMan.Init(angularRep-10, nAngularParameters-1);
365	// theBuff1.SetInterpolationManager(aMan); // Store interpolates f(costh)
366	// theBuff2.SetInterpolationManager(aMan); // Store interpolates f(costh)
367	for(i=0; i<nAngularParameters; i++) // i=1 ist wichtig!
368	{
369	theBuff1.SetX(i, theAngular[it-1].GetValue(i));
370	theBuff1.SetY(i, theAngular[it-1].GetValue(i+1));
371	theBuff2.SetX(i, theAngular[it].GetValue(i));
372	theBuff2.SetY(i, theAngular[it].GetValue(i+1));
373	i++;
374	}
375	G4NeutronHPVector theStore;
376	theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
377	x1 = y1;
378	x2 = y2;
379	G4double x, y;
380	//for(i=0;i<theBuff1.GetVectorLength(); i++);
381	for(i=0;i<theBuff1.GetVectorLength(); i++)
382	{
383	x = theBuff1.GetX(i); // costh binning identical
384	y1 = theBuff1.GetY(i);
385	y2 = theBuff2.GetY(i);
386	y = theInt.Interpolate(theManager.GetScheme(it),
387	fsEnergy, theAngular[it-1].GetLabel(),
388	theAngular[it].GetLabel(), y1, y2);
389	theStore.SetX(i, x);
390	theStore.SetY(i, y);
391	}
392	cosTh = theStore.Sample();
393	}
394	delete [] running;
395	}
396	else
397	{
398	throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPContAngularPar::Sample: Unknown angular representation");
399	}
400	result->SetKineticEnergy(fsEnergy);
401	G4double phi = twopi*G4UniformRand();
402	G4double theta = std::acos(cosTh);
403	G4double sinth = std::sin(theta);
404	G4double mtot = result->GetTotalMomentum();
405	G4ThreeVector tempVector(mtotsinthstd::cos(phi), mtotsinthstd::sin(phi), mtot*std::cos(theta) );
406	result->SetMomentum(tempVector);
407	// return the result.
408	return result;
409	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: