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