source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPElementData.cc@ 962

Last change on this file since 962 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 9.4 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// 02-08-06 Modified Harmonise to reslove cross section trouble at high-end. T. KOI
31//
32#include "G4NeutronHPElementData.hh"
33
34 G4NeutronHPElementData::G4NeutronHPElementData()
35 {
36 precision = 0.02;
37 theFissionData = new G4NeutronHPVector;
38 theCaptureData = new G4NeutronHPVector;
39 theElasticData = new G4NeutronHPVector;
40 theInelasticData = new G4NeutronHPVector;
41 theIsotopeWiseData = 0;
42 }
43
44 G4NeutronHPElementData::~G4NeutronHPElementData()
45 {
46 delete theFissionData;
47 delete theCaptureData;
48 delete theElasticData;
49 delete theInelasticData;
50 delete [] theIsotopeWiseData;
51 }
52
53 void G4NeutronHPElementData::Init(G4Element * theElement)
54 {
55 G4int count = theElement->GetNumberOfIsotopes();
56 if(count == 0) count +=
57 theStableOnes.GetNumberOfIsotopes(static_cast<G4int>(theElement->GetZ()));
58 theIsotopeWiseData = new G4NeutronHPIsoData[count];
59 // filename = ein data-set je isotope.
60 count = 0;
61 G4int nIso = theElement->GetNumberOfIsotopes();
62 G4int Z = static_cast<G4int> (theElement->GetZ());
63 G4int i1;
64 if(nIso!=0)
65 {
66 for (i1=0; i1<nIso; i1++)
67 {
68// G4cout <<" Init: normal case"<<G4endl;
69 G4int A = theElement->GetIsotope(i1)->GetN();
70 G4double frac = theElement->GetRelativeAbundanceVector()[i1]/perCent;
71 UpdateData(A, Z, count++, frac);
72 }
73 }else{
74// G4cout <<" Init: theStableOnes case: Z="<<Z<<G4endl;
75 G4int first = theStableOnes.GetFirstIsotope(Z);
76// G4cout <<"first="<<first<<" "<<theStableOnes.GetNumberOfIsotopes(theElement->GetZ())<<G4endl;
77 for(G4int i1=0;
78 i1<theStableOnes.GetNumberOfIsotopes(static_cast<G4int>(theElement->GetZ()) );
79 i1++)
80 {
81// G4cout <<" Init: theStableOnes in the loop"<<G4endl;
82 G4int A = theStableOnes.GetIsotopeNucleonCount(first+i1);
83 G4double frac = theStableOnes.GetAbundance(first+i1);
84// G4cout <<" Init: theStableOnes in the loop: "<<A<<G4endl;
85 UpdateData(A, Z, count++, frac);
86 }
87 }
88 theElasticData->ThinOut(precision);
89 theInelasticData->ThinOut(precision);
90 theCaptureData->ThinOut(precision);
91 theFissionData->ThinOut(precision);
92 }
93
94 void G4NeutronHPElementData::UpdateData(G4int A, G4int Z, G4int index, G4double abundance)
95 {
96 //Reads in the Data, using G4NeutronHPIsoData[], and its Init
97// G4cout << "entered: ElementWiseData::UpdateData"<<G4endl;
98 theIsotopeWiseData[index].Init(A, Z, abundance);
99// G4cout << "ElementWiseData::UpdateData Init finished"<<G4endl;
100
101 theBuffer = theIsotopeWiseData[index].MakeElasticData();
102// G4cout << "ElementWiseData::UpdateData MakeElasticData finished: "
103// <<theBuffer->GetVectorLength()<<G4endl;
104 Harmonise(theElasticData, theBuffer);
105// G4cout << "ElementWiseData::UpdateData Harmonise finished: "
106// <<theElasticData->GetVectorLength()<<G4endl;
107 delete theBuffer;
108
109 theBuffer = theIsotopeWiseData[index].MakeInelasticData();
110// G4cout << "ElementWiseData::UpdateData MakeInelasticData finished: "
111// <<theBuffer->GetVectorLength()<<G4endl;
112 Harmonise(theInelasticData, theBuffer);
113// G4cout << "ElementWiseData::UpdateData Harmonise finished: "
114// <<theInelasticData->GetVectorLength()<<G4endl;
115 delete theBuffer;
116
117 theBuffer = theIsotopeWiseData[index].MakeCaptureData();
118// G4cout << "ElementWiseData::UpdateData MakeCaptureData finished: "
119// <<theBuffer->GetVectorLength()<<G4endl;
120 Harmonise(theCaptureData, theBuffer);
121// G4cout << "ElementWiseData::UpdateData Harmonise finished: "
122// <<theCaptureData->GetVectorLength()<<G4endl;
123 delete theBuffer;
124
125 theBuffer = theIsotopeWiseData[index].MakeFissionData();
126// G4cout << "ElementWiseData::UpdateData MakeFissionData finished: "
127// <<theBuffer->GetVectorLength()<<G4endl;
128 Harmonise(theFissionData, theBuffer);
129// G4cout << "ElementWiseData::UpdateData Harmonise finished: "
130// <<theFissionData->GetVectorLength()<<G4endl;
131 delete theBuffer;
132
133// G4cout << "ElementWiseData::UpdateData finished"<endl;
134 }
135
136 void G4NeutronHPElementData::Harmonise(G4NeutronHPVector *& theStore, G4NeutronHPVector * theNew)
137 {
138 if(theNew == 0) { return; }
139 G4int s = 0, n=0, m=0;
140 G4NeutronHPVector * theMerge = new G4NeutronHPVector(theStore->GetVectorLength());
141// G4cout << "Harmonise 1: "<<theStore->GetEnergy(s)<<" "<<theNew->GetEnergy(0)<<G4endl;
142 while ( theStore->GetEnergy(s)<theNew->GetEnergy(0)&&s<theStore->GetVectorLength() )
143 {
144 theMerge->SetData(m++, theStore->GetEnergy(s), theStore->GetXsec(s));
145 s++;
146 }
147 G4NeutronHPVector *active = theStore;
148 G4NeutronHPVector * passive = theNew;
149 G4NeutronHPVector * tmp;
150 G4int a = s, p = n, t;
151// G4cout << "Harmonise 2: "<<active->GetVectorLength()<<" "<<passive->GetVectorLength()<<G4endl;
152 while (a<active->GetVectorLength()&&p<passive->GetVectorLength())
153 {
154 if(active->GetEnergy(a) <= passive->GetEnergy(p))
155 {
156 theMerge->SetData(m, active->GetEnergy(a), active->GetXsec(a));
157 G4double x = theMerge->GetEnergy(m);
158 G4double y = std::max(0., passive->GetXsec(x));
159 theMerge->SetData(m, x, theMerge->GetXsec(m)+y);
160 m++;
161 a++;
162 } else {
163// G4cout << "swapping in Harmonise"<<G4endl;
164 tmp = active; t=a;
165 active = passive; a=p;
166 passive = tmp; p=t;
167 }
168 }
169// G4cout << "Harmonise 3: "<< a <<" "<<active->GetVectorLength()<<" "<<m<<G4endl;
170 while (a!=active->GetVectorLength())
171 {
172 theMerge->SetData(m++, active->GetEnergy(a), active->GetXsec(a));
173 a++;
174 }
175// G4cout << "Harmonise 4: "<< p <<" "<<passive->GetVectorLength()<<" "<<m<<G4endl;
176 while (p!=passive->GetVectorLength())
177 {
178 // Modified by T. KOI
179 //theMerge->SetData(m++, passive->GetEnergy(p), passive->GetXsec(p));
180 G4double x = passive->GetEnergy(p);
181 G4double y = std::max(0., active->GetXsec(x));
182 theMerge->SetData(m++, x, passive->GetXsec(p)+y);
183 p++;
184 }
185// G4cout <<"Harmonise 5: "<< theMerge->GetVectorLength() << " " << m << G4endl;
186 delete theStore;
187 theStore = theMerge;
188// G4cout <<"Harmonise 6: "<< theStore->GetVectorLength() << " " << m << G4endl;
189 }
190
191 G4NeutronHPVector * G4NeutronHPElementData::MakePhysicsVector(G4Element * theElement,
192 G4ParticleDefinition * theP,
193 G4NeutronHPFissionData* theSet)
194 {
195 if(theP != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
196 Init ( theElement );
197 return GetData(theSet);
198 }
199 G4NeutronHPVector * G4NeutronHPElementData::MakePhysicsVector(G4Element * theElement,
200 G4ParticleDefinition * theP,
201 G4NeutronHPCaptureData * theSet)
202 {
203 if(theP != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
204 Init ( theElement );
205 return GetData(theSet);
206 }
207 G4NeutronHPVector * G4NeutronHPElementData::MakePhysicsVector(G4Element * theElement,
208 G4ParticleDefinition * theP,
209 G4NeutronHPElasticData * theSet)
210 {
211 if(theP != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
212 Init ( theElement );
213 return GetData(theSet);
214 }
215 G4NeutronHPVector * G4NeutronHPElementData::MakePhysicsVector(G4Element * theElement,
216 G4ParticleDefinition * theP,
217 G4NeutronHPInelasticData * theSet)
218 {
219 if(theP != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
220 Init ( theElement );
221 return GetData(theSet);
222 }
Note: See TracBrowser for help on using the repository browser.