Context Navigation

G4NeutronHPContAngularPar.cc @ 968

Last change on this file since 968 was 962, checked in by garnier, 15 years ago
update processes
File size: 18.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	// 080612 Fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #1
35	// 080714 Limiting the sum of energy of secondary particles by T. Koi
36	// 080801 Fix div0 error wiht G4FPE and memory leak by T. Koi
37	// 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
38	//
39
40	#include "G4NeutronHPContAngularPar.hh"
41	#include "G4NeutronHPLegendreStore.hh"
42	#include "G4Gamma.hh"
43	#include "G4Electron.hh"
44	#include "G4Positron.hh"
45	#include "G4Neutron.hh"
46	#include "G4Proton.hh"
47	#include "G4Deuteron.hh"
48	#include "G4Triton.hh"
49	#include "G4He3.hh"
50	#include "G4Alpha.hh"
51	#include "G4NeutronHPVector.hh"
52	#include "G4NucleiProperties.hh"
53	#include "G4NeutronHPKallbachMannSyst.hh"
54	#include "G4ParticleTable.hh"
55
56	void G4NeutronHPContAngularPar::Init(std::ifstream & aDataFile)
57	{
58	aDataFile >> theEnergy >> nEnergies >> nDiscreteEnergies >> nAngularParameters;
59	theEnergy *= eV;
60	theAngular = new G4NeutronHPList [nEnergies];
61	for(G4int i=0; i<nEnergies; i++)
62	{
63	G4double sEnergy;
64	aDataFile >> sEnergy;
65	sEnergy*=eV;
66	theAngular[i].SetLabel(sEnergy);
67	theAngular[i].Init(aDataFile, nAngularParameters, 1.);
68	}
69	}
70
71	G4ReactionProduct *
72	G4NeutronHPContAngularPar::Sample(G4double anEnergy, G4double massCode, G4double /targetMass/,
73	G4int angularRep, G4int interpolE )
74	{
75	G4ReactionProduct * result = new G4ReactionProduct;
76	G4int Z = static_cast<G4int>(massCode/1000);
77	G4int A = static_cast<G4int>(massCode-1000*Z);
78	if(massCode==0)
79	{
80	result->SetDefinition(G4Gamma::Gamma());
81	}
82	else if(A==0)
83	{
84	result->SetDefinition(G4Electron::Electron());
85	if(Z==1) result->SetDefinition(G4Positron::Positron());
86	}
87	else if(A==1)
88	{
89	result->SetDefinition(G4Neutron::Neutron());
90	if(Z==1) result->SetDefinition(G4Proton::Proton());
91	}
92	else if(A==2)
93	{
94	result->SetDefinition(G4Deuteron::Deuteron());
95	}
96	else if(A==3)
97	{
98	result->SetDefinition(G4Triton::Triton());
99	if(Z==2) result->SetDefinition(G4He3::He3());
100	}
101	else if(A==4)
102	{
103	result->SetDefinition(G4Alpha::Alpha());
104	if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPContAngularPar: Unknown ion case 1");
105	}
106	else
107	{
108	result->SetDefinition(G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z));
109	}
110	G4int i(0);
111	G4int it(0);
112	G4double fsEnergy(0);
113	G4double cosTh(0);
114	if(angularRep==1)
115	{
116	// 080612 Fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #1
117	if (interpolE == 2)
118	{
119
120	//TK080711
121	if ( fresh == true )
122	{
123	remaining_energy = theAngular[0].GetLabel();
124	fresh = false;
125	}
126	//TK080711
127
128	G4double random = G4UniformRand();
129	G4double * running = new G4double[nEnergies+1];
130	running[0]=0;
131
132	for(i=1; i<nEnergies+1; i++)
133	{
134	//TK080711
135	if ( remaining_energy >= theAngular[ i-1 ].GetLabel() )
136	running[i] = running[i-1] + theAngular[i-1].GetValue(0);
137	else
138	running[i] = running[i-1];
139	//TK080711
140	}
141
142	//080730
143	if ( running[ nEnergies ] != 0 )
144	{
145
146	for ( i = 1 ; i < nEnergies+1 ; i++ )
147	{
148	it = i-1;
149	if ( random > running[ i-1 ]/running[ nEnergies ] && random <= running[ i ] / running[ nEnergies ] ) break;
150	}
151	fsEnergy = theAngular[ it ].GetLabel();
152
153	}
154
155	//TK080711
156	if ( i == nEnergies+1 \|\| running[ nEnergies ] == 0 ) fsEnergy = remaining_energy;
157	//TK080711 //080730
158
159	G4NeutronHPLegendreStore theStore(1);
160	theStore.Init(0,fsEnergy,nAngularParameters);
161	for(i=0;i<nAngularParameters;i++)
162	{
163	theStore.SetCoeff(0,i,theAngular[it].GetValue(i));
164	}
165	// use it to sample.
166	cosTh = theStore.SampleMax(fsEnergy);
167
168	//TK080711
169	remaining_energy -= fsEnergy;
170	//TK080711
171
172	//080801b
173	delete[] running;
174	//080801b
175	}
176	else
177	{
178
179	//080714
180	if ( fresh == true )
181	{
182	remaining_energy = theAngular[ nEnergies-1 ].GetLabel();
183	fresh = false;
184	}
185	//080714
186
187
188	G4double random = G4UniformRand();
189	G4double * running = new G4double[nEnergies];
190	running[0]=0;
191	G4double weighted = 0;
192	for(i=1; i<nEnergies; i++)
193	{
194	/*
195	if(i!=0)
196	{
197	running[i]=running[i-1];
198	}
199	running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
200	theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
201	theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
202	weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
203	theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
204	theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
205	*/
206
207	running[i]=running[i-1];
208	if ( remaining_energy >= theAngular[i].GetLabel() )
209	{
210	running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
211	theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
212	theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
213	weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
214	theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
215	theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
216	}
217	}
218	// cash the mean energy in this distribution
219	//080409 TKDB
220	if ( nEnergies == 1 \|\| running[nEnergies-1] == 0 )
221	currentMeanEnergy = 0.0;
222	else
223	{
224	currentMeanEnergy = weighted/running[nEnergies-1];
225	}
226
227	//080409 TKDB
228	if ( nEnergies == 1 ) it = 0;
229
230	//080729
231	if ( running[nEnergies-1] != 0 )
232	{
233	for ( i = 1 ; i < nEnergies ; i++ )
234	{
235	it = i;
236	if ( random < running [ i ] / running [ nEnergies-1 ] ) break;
237	}
238	}
239
240	//080714
241	if ( running [ nEnergies-1 ] == 0 ) it = 0;
242	//080714
243
244	if(it<nDiscreteEnergies\|\|it==0)
245	{
246	if(it == 0)
247	{
248	fsEnergy = theAngular[it].GetLabel();
249	G4NeutronHPLegendreStore theStore(1);
250	theStore.Init(0,fsEnergy,nAngularParameters);
251	for(i=0;i<nAngularParameters;i++)
252	{
253	theStore.SetCoeff(0,i,theAngular[it].GetValue(i));
254	}
255	// use it to sample.
256	cosTh = theStore.SampleMax(fsEnergy);
257	}
258	else
259	{
260	G4double e1, e2;
261	e1 = theAngular[it-1].GetLabel();
262	e2 = theAngular[it].GetLabel();
263	fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
264	random,
265	running[it-1]/running[nEnergies-1],
266	running[it]/running[nEnergies-1],
267	e1, e2);
268	// fill a Legendrestore
269	G4NeutronHPLegendreStore theStore(2);
270	theStore.Init(0,e1,nAngularParameters);
271	theStore.Init(1,e2,nAngularParameters);
272	for(i=0;i<nAngularParameters;i++)
273	{
274	theStore.SetCoeff(0,i,theAngular[it-1].GetValue(i));
275	theStore.SetCoeff(1,i,theAngular[it].GetValue(i));
276	}
277	// use it to sample.
278	theStore.SetManager(theManager);
279	cosTh = theStore.SampleMax(fsEnergy);
280	}
281	}
282	else // continuum contribution
283	{
284	G4double x1 = running[it-1]/running[nEnergies-1];
285	G4double x2 = running[it]/running[nEnergies-1];
286	G4double y1 = theAngular[it-1].GetLabel();
287	G4double y2 = theAngular[it].GetLabel();
288	fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
289	random,x1,x2,y1,y2);
290	G4NeutronHPLegendreStore theStore(2);
291	theStore.Init(0,y1,nAngularParameters);
292	theStore.Init(1,y2,nAngularParameters);
293	theStore.SetManager(theManager);
294	for(i=0;i<nAngularParameters;i++)
295	{
296	theStore.SetCoeff(0,i,theAngular[it-1].GetValue(i));
297	theStore.SetCoeff(1,i,theAngular[it].GetValue(i));
298	}
299	// use it to sample.
300	cosTh = theStore.SampleMax(fsEnergy);
301	}
302	delete [] running;
303
304	//080714
305	remaining_energy -= fsEnergy;
306	//080714
307
308	}
309
310	}
311	else if(angularRep==2)
312	{
313	// first get the energy (already the right for this incoming energy)
314	G4int i;
315	G4double * running = new G4double[nEnergies];
316	running[0]=0;
317	G4double weighted = 0;
318	for(i=1; i<nEnergies; i++)
319	{
320	if(i!=0) running[i]=running[i-1];
321	running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
322	theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
323	theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
324	weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
325	theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
326	theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
327	}
328	// cash the mean energy in this distribution
329	//080409 TKDB
330	//currentMeanEnergy = weighted/running[nEnergies-1];
331	if ( nEnergies == 1 )
332	currentMeanEnergy = 0.0;
333	else
334	currentMeanEnergy = weighted/running[nEnergies-1];
335
336	G4int it(0);
337	G4double randkal = G4UniformRand();
338	//080409 TKDB
339	//for(i=0; i<nEnergies; i++)
340	for(i=1; i<nEnergies; i++)
341	{
342	it = i;
343	if(randkal<running[i]/running[nEnergies-1]) break;
344	}
345
346	// interpolate the secondary energy.
347	G4double x, x1,x2,y1,y2;
348	if(it==0) it=1;
349	x = randkal*running[nEnergies-1];
350	x1 = running[it-1];
351	x2 = running[it];
352	G4double compoundFraction;
353	// interpolate energy
354	y1 = theAngular[it-1].GetLabel();
355	y2 = theAngular[it].GetLabel();
356	fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it-1),
357	x, x1,x2,y1,y2);
358	// for theta interpolate the compoundFractions
359	G4double cLow = theAngular[it-1].GetValue(1);
360	G4double cHigh = theAngular[it].GetValue(1);
361	compoundFraction = theInt.Interpolate(theManager.GetScheme(it),
362	fsEnergy, y1, y2, cLow,cHigh);
363	delete [] running;
364
365	// get cosTh
366	G4double incidentEnergy = anEnergy;
367	G4double incidentMass = G4Neutron::Neutron()->GetPDGMass();
368	G4double productEnergy = fsEnergy;
369	G4double productMass = result->GetMass();
370	G4int targetZ = G4int(theTargetCode/1000);
371	G4int targetA = G4int(theTargetCode-1000*targetZ);
372	// To correspond to natural composition (-nat-) data files.
373	if ( targetA == 0 )
374	targetA = int ( theTarget->GetMass()/amu_c2 + 0.5 );
375	G4double targetMass = theTarget->GetMass();
376	G4int residualA = targetA+1-A;
377	G4int residualZ = targetZ-Z;
378	G4double residualMass = residualZ*G4Proton::Proton()->GetPDGMass();
379	residualMass +=(residualA-residualZ)*G4Neutron::Neutron()->GetPDGMass();
380	residualMass -= G4NucleiProperties::GetBindingEnergy( residualA , residualZ );
381	G4NeutronHPKallbachMannSyst theKallbach(compoundFraction,
382	incidentEnergy, incidentMass,
383	productEnergy, productMass,
384	residualMass, residualA, residualZ,
385	targetMass, targetA, targetZ);
386	cosTh = theKallbach.Sample(anEnergy);
387	}
388	else if(angularRep>10&&angularRep<16)
389	{
390	G4double random = G4UniformRand();
391	G4double * running = new G4double[nEnergies];
392	running[0]=0;
393	G4double weighted = 0;
394	for(i=1; i<nEnergies; i++)
395	{
396	if(i!=0) running[i]=running[i-1];
397	running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
398	theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
399	theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
400	weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
401	theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
402	theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
403	}
404	// cash the mean energy in this distribution
405	//currentMeanEnergy = weighted/running[nEnergies-1];
406	if ( nEnergies == 1 )
407	currentMeanEnergy = 0.0;
408	else
409	currentMeanEnergy = weighted/running[nEnergies-1];
410
411	//080409 TKDB
412	if ( nEnergies == 1 ) it = 0;
413	//for(i=0; i<nEnergies; i++)
414	for(i=1; i<nEnergies; i++)
415	{
416	it = i;
417	if(random<running[i]/running[nEnergies-1]) break;
418	}
419	if(it<nDiscreteEnergies\|\|it==0)
420	{
421	if(it==0)
422	{
423	fsEnergy = theAngular[0].GetLabel();
424	G4NeutronHPVector theStore;
425	G4int aCounter = 0;
426	for(G4int i=1; i<nAngularParameters; i+=2)
427	{
428	theStore.SetX(aCounter, theAngular[0].GetValue(i));
429	theStore.SetY(aCounter, theAngular[0].GetValue(i+1));
430	aCounter++;
431	}
432	G4InterpolationManager aMan;
433	aMan.Init(angularRep-10, nAngularParameters-1);
434	theStore.SetInterpolationManager(aMan);
435	cosTh = theStore.Sample();
436	}
437	else
438	{
439	fsEnergy = theAngular[it].GetLabel();
440	G4NeutronHPVector theStore;
441	G4InterpolationManager aMan;
442	aMan.Init(angularRep-10, nAngularParameters-1);
443	theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
444	G4InterpolationScheme currentScheme = theManager.GetInverseScheme(it);
445	G4int aCounter = 0;
446	for(G4int i=1; i<nAngularParameters; i+=2)
447	{
448	theStore.SetX(aCounter, theAngular[it].GetValue(i));
449	theStore.SetY(aCounter, theInt.Interpolate(currentScheme,
450	random,
451	running[it-1]/running[nEnergies-1],
452	running[it]/running[nEnergies-1],
453	theAngular[it-1].GetValue(i+1),
454	theAngular[it].GetValue(i+1)));
455	aCounter++;
456	}
457	cosTh = theStore.Sample();
458	}
459	}
460	else
461	{
462	G4double x1 = running[it-1]/running[nEnergies-1];
463	G4double x2 = running[it]/running[nEnergies-1];
464	G4double y1 = theAngular[it-1].GetLabel();
465	G4double y2 = theAngular[it].GetLabel();
466	fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
467	random,x1,x2,y1,y2);
468	G4NeutronHPVector theBuff1;
469	G4NeutronHPVector theBuff2;
470	G4InterpolationManager aMan;
471	aMan.Init(angularRep-10, nAngularParameters-1);
472	// theBuff1.SetInterpolationManager(aMan); // Store interpolates f(costh)
473	// theBuff2.SetInterpolationManager(aMan); // Store interpolates f(costh)
474	for(i=0; i<nAngularParameters; i++) // i=1 ist wichtig!
475	{
476	theBuff1.SetX(i, theAngular[it-1].GetValue(i));
477	theBuff1.SetY(i, theAngular[it-1].GetValue(i+1));
478	theBuff2.SetX(i, theAngular[it].GetValue(i));
479	theBuff2.SetY(i, theAngular[it].GetValue(i+1));
480	i++;
481	}
482	G4NeutronHPVector theStore;
483	theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
484	x1 = y1;
485	x2 = y2;
486	G4double x, y;
487	//for(i=0;i<theBuff1.GetVectorLength(); i++);
488	for(i=0;i<theBuff1.GetVectorLength(); i++)
489	{
490	x = theBuff1.GetX(i); // costh binning identical
491	y1 = theBuff1.GetY(i);
492	y2 = theBuff2.GetY(i);
493	y = theInt.Interpolate(theManager.GetScheme(it),
494	fsEnergy, theAngular[it-1].GetLabel(),
495	theAngular[it].GetLabel(), y1, y2);
496	theStore.SetX(i, x);
497	theStore.SetY(i, y);
498	}
499	cosTh = theStore.Sample();
500	}
501	delete [] running;
502	}
503	else
504	{
505	throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPContAngularPar::Sample: Unknown angular representation");
506	}
507	result->SetKineticEnergy(fsEnergy);
508	G4double phi = twopi*G4UniformRand();
509	G4double theta = std::acos(cosTh);
510	G4double sinth = std::sin(theta);
511	G4double mtot = result->GetTotalMomentum();
512	G4ThreeVector tempVector(mtotsinthstd::cos(phi), mtotsinthstd::sin(phi), mtot*std::cos(theta) );
513	result->SetMomentum(tempVector);
514	// return the result.
515	return result;
516	}

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

Context Navigation

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

Download in other formats: