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

Last change on this file since 1340 was 819, checked in by garnier, 16 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.