source: trunk/source/processes/hadronic/cross_sections/src/G4CrossSectionDataStore.cc @ 846

Last change on this file since 846 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 10.8 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//
27
28#include "G4CrossSectionDataStore.hh"
29#include "G4HadronicException.hh"
30#include "G4StableIsotopes.hh"
31#include "G4HadTmpUtil.hh"
32#include "Randomize.hh"
33
34
35G4double
36G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* aParticle,
37                                         const G4Element* anElement,
38                                         G4double aTemperature)
39{
40  if (NDataSetList == 0) 
41  {
42    throw G4HadronicException(__FILE__, __LINE__, 
43     "G4CrossSectionDataStore: no data sets registered");
44    return DBL_MIN;
45  }
46  for (G4int i = NDataSetList-1; i >= 0; i--) {
47    if (DataSetList[i]->IsApplicable(aParticle, anElement))
48      return DataSetList[i]->GetCrossSection(aParticle,anElement,aTemperature);
49  }
50  throw G4HadronicException(__FILE__, __LINE__, 
51                      "G4CrossSectionDataStore: no applicable data set found "
52                      "for particle/element");
53  return DBL_MIN;
54}
55
56G4VCrossSectionDataSet*
57G4CrossSectionDataStore::whichDataSetInCharge(const G4DynamicParticle* aParticle,
58                                              const G4Element* anElement)
59{
60  if (NDataSetList == 0) 
61  {
62    throw G4HadronicException(__FILE__, __LINE__, 
63     "G4CrossSectionDataStore: no data sets registered");
64    return 0;
65  }
66  for (G4int i = NDataSetList-1; i >= 0; i--) {
67    if (DataSetList[i]->IsApplicable(aParticle, anElement) )
68    {
69      return DataSetList[i];
70    }
71  }
72  throw G4HadronicException(__FILE__, __LINE__, 
73                      "G4CrossSectionDataStore: no applicable data set found "
74                      "for particle/element");
75  return 0;
76}
77
78
79G4double
80G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* aParticle,
81                                         const G4Isotope* anIsotope,
82                                         G4double aTemperature)
83{
84  if (NDataSetList == 0) 
85  {
86    throw G4HadronicException(__FILE__, __LINE__, 
87     "G4CrossSectionDataStore: no data sets registered");
88    return DBL_MIN;
89  }
90  for (G4int i = NDataSetList-1; i >= 0; i--) {
91    if (DataSetList[i]->IsZAApplicable(aParticle, anIsotope->GetZ(), anIsotope->GetN()))
92      return DataSetList[i]->GetIsoCrossSection(aParticle,anIsotope,aTemperature);
93  }
94  throw G4HadronicException(__FILE__, __LINE__, 
95                      "G4CrossSectionDataStore: no applicable data set found "
96                      "for particle/element");
97  return DBL_MIN;
98}
99
100
101G4double
102G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* aParticle,
103                                         G4double Z, G4double A,
104                                         G4double aTemperature)
105{
106  if (NDataSetList == 0) 
107  {
108    throw G4HadronicException(__FILE__, __LINE__, 
109     "G4CrossSectionDataStore: no data sets registered");
110    return DBL_MIN;
111  }
112  for (G4int i = NDataSetList-1; i >= 0; i--) {
113      if (DataSetList[i]->IsZAApplicable(aParticle, Z, A))
114      return DataSetList[i]->GetIsoZACrossSection(aParticle,Z,A,aTemperature);
115  }
116  throw G4HadronicException(__FILE__, __LINE__, 
117                      "G4CrossSectionDataStore: no applicable data set found "
118                      "for particle/element");
119  return DBL_MIN;
120}
121
122
123G4double
124G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* aParticle,
125                                         const G4Material* aMaterial)
126{
127  G4double sigma(0);
128  G4double aTemp = aMaterial->GetTemperature();
129  G4int nElements = aMaterial->GetNumberOfElements();
130  const G4double* theAtomsPerVolumeVector = aMaterial->GetVecNbOfAtomsPerVolume();
131  G4double xSection(0);
132
133  for(G4int i=0; i<nElements; i++) {
134    xSection = GetCrossSection( aParticle, (*aMaterial->GetElementVector())[i], aTemp);
135    sigma += theAtomsPerVolumeVector[i] * xSection;
136  }
137
138  return sigma;
139}
140
141
142std::pair<G4double, G4double> 
143G4CrossSectionDataStore::SelectRandomIsotope(const G4DynamicParticle* particle,
144                                             const G4Material* aMaterial)
145{
146  static G4StableIsotopes theDefaultIsotopes;  // natural abundances and
147                                               // stable isotopes
148  G4double aTemp = aMaterial->GetTemperature();
149  G4int nElements = aMaterial->GetNumberOfElements();
150  const G4ElementVector* theElementVector = aMaterial->GetElementVector();
151  const G4double* theAtomsPerVolumeVector = aMaterial->GetVecNbOfAtomsPerVolume();
152  std::vector<std::vector<G4double> > awicsPerElement;
153  std::vector<std::vector<G4double> > AvaluesPerElement;
154  G4Element* anElement;
155
156  // Collect abundance weighted cross sections and A values for each isotope
157  // in each element
158 
159  for (G4int i = 0; i < nElements; i++) {
160    anElement = (*theElementVector)[i];
161    G4int nIsoPerElement = anElement->GetNumberOfIsotopes();
162    std::vector<G4double> isoholder;
163    std::vector<G4double> aholder;
164    G4double iso_xs = DBL_MIN;
165
166    if (nIsoPerElement) { // user-defined isotope abundances
167      G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
168      G4double* abundVector = anElement->GetRelativeAbundanceVector();
169      G4VCrossSectionDataSet* inCharge = whichDataSetInCharge(particle, anElement);
170      G4bool elementXS = false;
171      for (G4int j = 0; j < nIsoPerElement; j++) {
172        if (inCharge->IsZAApplicable(particle, (*isoVector)[j]->GetZ(), 
173                                               (*isoVector)[j]->GetN() ) ) {
174          iso_xs = inCharge->GetIsoCrossSection(particle, (*isoVector)[j], aTemp);
175        } else if (elementXS == false) {
176          iso_xs = inCharge->GetCrossSection(particle, anElement, aTemp);
177          elementXS = true;
178        }
179
180        isoholder.push_back(abundVector[j]*iso_xs);
181        aholder.push_back(G4double((*isoVector)[j]->GetN()));
182      }
183
184    } else { // natural abundances
185      G4int ZZ = G4lrint(anElement->GetZ());
186      nIsoPerElement = theDefaultIsotopes.GetNumberOfIsotopes(ZZ);
187      G4int index = theDefaultIsotopes.GetFirstIsotope(ZZ);
188      G4double AA;
189      G4double abundance;
190
191      G4VCrossSectionDataSet* inCharge = whichDataSetInCharge(particle, anElement);
192      G4bool elementXS = false;
193
194      for (G4int j = 0; j < nIsoPerElement; j++) {
195        AA = G4double(theDefaultIsotopes.GetIsotopeNucleonCount(index+j));
196        aholder.push_back(AA);
197        if (inCharge->IsZAApplicable(particle, G4double(ZZ), AA) ) {
198          iso_xs = inCharge->GetIsoZACrossSection(particle, G4double(ZZ), AA, aTemp);
199        } else if (elementXS == false) {
200          iso_xs = inCharge->GetCrossSection(particle, anElement, aTemp);
201          elementXS = true;
202        }
203        abundance = theDefaultIsotopes.GetAbundance(index+j)/100.0;
204        isoholder.push_back(abundance*iso_xs);
205      }
206    }
207
208    awicsPerElement.push_back(isoholder);
209    AvaluesPerElement.push_back(aholder);
210  }
211
212  // Calculate running sums for isotope selection
213
214  G4double crossSectionTotal = 0;
215  G4double xSectionPerElement;
216  std::vector<G4double> runningSum;
217
218  for (G4int i=0; i < nElements; i++) {
219    xSectionPerElement = 0;
220    for (G4int j=0; j < G4int(awicsPerElement[i].size()); j++)
221                     xSectionPerElement += awicsPerElement[i][j];
222    runningSum.push_back(theAtomsPerVolumeVector[i]*xSectionPerElement);
223    crossSectionTotal += runningSum[i];
224  }
225
226  // Compare random number to running sum over element xc to choose Z
227
228  // Initialize Z and A to first element and first isotope in case
229  // cross section is zero
230   
231  G4double ZZ = (*theElementVector)[0]->GetZ();
232  G4double AA = AvaluesPerElement[0][0];
233  if (crossSectionTotal != 0.) {
234    G4double random = G4UniformRand();
235    for(G4int i=0; i < nElements; i++) {
236      if(i!=0) runningSum[i] += runningSum[i-1];
237      if(random <= runningSum[i]/crossSectionTotal) {
238        ZZ = ((*theElementVector)[i])->GetZ();
239
240        // Compare random number to running sum over isotope xc to choose A
241
242        G4int nIso = awicsPerElement[i].size();
243        G4double* running = new G4double[nIso];
244        for (G4int j=0; j < nIso; j++) {
245          running[j] = awicsPerElement[i][j];
246          if(j!=0) running[j] += running[j-1];
247        }
248
249        G4double trial = G4UniformRand(); 
250        for (G4int j=0; j < nIso; j++) {
251          AA = AvaluesPerElement[i][j];
252          if (trial <= running[j]/running[nIso-1]) break;
253        }
254        delete [] running;
255        break;
256      }
257    }
258  }
259  return std::pair<G4double, G4double>(ZZ, AA);
260}
261
262
263void
264G4CrossSectionDataStore::AddDataSet(G4VCrossSectionDataSet* aDataSet)
265{
266   if (NDataSetList == NDataSetMax) {
267      G4cout << "WARNING: G4CrossSectionDataStore::AddDataSet: "<<G4endl;
268      G4cout << "         reached maximum number of data sets";
269      G4cout << "         data set not added !!!!!!!!!!!!!!!!";
270      return;
271   }
272   DataSetList[NDataSetList] = aDataSet;
273   NDataSetList++;
274}
275
276
277void
278G4CrossSectionDataStore::
279BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
280{
281   if (NDataSetList == 0) 
282   {
283     G4Exception("G4CrossSectionDataStore", "007", FatalException,
284                 "BuildPhysicsTable: no data sets registered");
285     return;
286   }
287   for (G4int i = NDataSetList-1; i >= 0; i--) {
288      DataSetList[i]->BuildPhysicsTable(aParticleType);
289   }
290}
291
292
293void
294G4CrossSectionDataStore::
295DumpPhysicsTable(const G4ParticleDefinition& aParticleType)
296{
297   if (NDataSetList == 0) {
298      G4cout << "WARNING - G4CrossSectionDataStore::DumpPhysicsTable: no data sets registered"<<G4endl;
299      return;
300   }
301   for (G4int i = NDataSetList-1; i >= 0; i--) {
302      DataSetList[i]->DumpPhysicsTable(aParticleType);
303   }
304}
Note: See TracBrowser for help on using the repository browser.