Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4CrossSectionDataStore.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 // $Id$
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4CrossSectionDataStore
34 //
35 // Modifications:
36 // 23.01.2009 V.Ivanchenko add destruction of data sets
37 // 29.04.2010 G.Folger modifictaions for integer A & Z
38 // 14.03.2011 V.Ivanchenko fixed DumpPhysicsTable
39 // 15.08.2011 G.Folger, V.Ivanchenko, T.Koi, D.Wright redesign the class
40 //
41 
42 #include "G4CrossSectionDataStore.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4HadronicException.hh"
45 #include "G4HadTmpUtil.hh"
46 #include "Randomize.hh"
47 #include "G4Nucleus.hh"
48 
49 #include "G4DynamicParticle.hh"
50 #include "G4Isotope.hh"
51 #include "G4Element.hh"
52 #include "G4Material.hh"
53 #include "G4NistManager.hh"
54 #include <iostream>
55 
56 
57 G4CrossSectionDataStore::G4CrossSectionDataStore() :
58  nDataSetList(0), verboseLevel(0)
59 {
60  nist = G4NistManager::Instance();
61  currentMaterial = elmMaterial = 0;
62  currentElement = 0; //ALB 14-Aug-2012 Coverity fix.
63  matParticle = elmParticle = 0;
64  matKinEnergy = elmKinEnergy = matCrossSection = elmCrossSection = 0.0;
65 }
66 
67 G4CrossSectionDataStore::~G4CrossSectionDataStore()
68 {}
69 
70 G4double
71 G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* part,
72  const G4Material* mat)
73 {
74  if(mat == currentMaterial && part->GetDefinition() == matParticle
75  && part->GetKineticEnergy() == matKinEnergy)
76  { return matCrossSection; }
77 
78  currentMaterial = mat;
79  matParticle = part->GetDefinition();
80  matKinEnergy = part->GetKineticEnergy();
81  matCrossSection = 0;
82 
83  G4int nElements = mat->GetNumberOfElements();
84  const G4double* nAtomsPerVolume = mat->GetVecNbOfAtomsPerVolume();
85 
86  if(G4int(xsecelm.size()) < nElements) { xsecelm.resize(nElements); }
87 
88  for(G4int i=0; i<nElements; ++i) {
89  matCrossSection += nAtomsPerVolume[i] *
90  GetCrossSection(part, (*mat->GetElementVector())[i], mat);
91  xsecelm[i] = matCrossSection;
92  }
93  return matCrossSection;
94 }
95 
96 G4double
97 G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* part,
98  const G4Element* elm,
99  const G4Material* mat)
100 {
101  if(mat == elmMaterial && elm == currentElement &&
102  part->GetDefinition() == elmParticle &&
103  part->GetKineticEnergy() == elmKinEnergy)
104  { return elmCrossSection; }
105 
106  elmMaterial = mat;
107  currentElement = elm;
108  elmParticle = part->GetDefinition();
109  elmKinEnergy = part->GetKineticEnergy();
110  elmCrossSection = 0.0;
111 
112  G4int i = nDataSetList-1;
113  G4int Z = G4lrint(elm->GetZ());
114  if (dataSetList[i]->IsElementApplicable(part, Z, mat)) {
115 
116  // element wise cross section
117  elmCrossSection = dataSetList[i]->GetElementCrossSection(part, Z, mat);
118 
119  } else {
120  // isotope wise cross section
121  G4int nIso = elm->GetNumberOfIsotopes();
122  G4Isotope* iso = 0;
123 
124  if (0 >= nIso) {
125  G4cout << " Element " << elm->GetName() << " Z= " << Z
126  << " has no isotopes " << G4endl;
127  throw G4HadronicException(__FILE__, __LINE__,
128  " Isotope vector is not defined");
129  //ALB 14-Aug-2012 Coverity fix. return elmCrossSection;
130  }
131  // user-defined isotope abundances
132  G4IsotopeVector* isoVector = elm->GetIsotopeVector();
133  G4double* abundVector = elm->GetRelativeAbundanceVector();
134 
135  for (G4int j = 0; j<nIso; ++j) {
136  if(abundVector[j] > 0.0) {
137  iso = (*isoVector)[j];
138  elmCrossSection += abundVector[j]*
139  GetIsoCrossSection(part, Z, iso->GetN(), iso, elm, mat, i);
140  }
141  }
142  }
143  //G4cout << "xsec(barn)= " << xsec/barn << G4endl;
144  return elmCrossSection;
145 }
146 
147 G4double
148 G4CrossSectionDataStore::GetIsoCrossSection(const G4DynamicParticle* part,
149  G4int Z, G4int A,
150  const G4Isotope* iso,
151  const G4Element* elm,
152  const G4Material* mat,
153  G4int idx)
154 {
155  // this methods is called after the check that dataSetList[idx]
156  // depend on isotopes, so for this DataSet only isotopes are checked
157 
158  // isotope-wise cross section does exist
159  if(dataSetList[idx]->IsIsoApplicable(part, Z, A, elm, mat) ) {
160  return dataSetList[idx]->GetIsoCrossSection(part, Z, A, iso, elm, mat);
161 
162  } else {
163  // seach for other dataSet
164  for (G4int j = idx-1; j >= 0; --j) {
165  if (dataSetList[j]->IsElementApplicable(part, Z, mat)) {
166  return dataSetList[j]->GetElementCrossSection(part, Z, mat);
167  } else if (dataSetList[j]->IsIsoApplicable(part, Z, A, elm, mat)) {
168  return dataSetList[j]->GetIsoCrossSection(part, Z, A, iso, elm, mat);
169  }
170  }
171  }
172  G4cout << "G4CrossSectionDataStore::GetCrossSection ERROR: "
173  << " no isotope cross section found"
174  << G4endl;
175  G4cout << " for " << part->GetDefinition()->GetParticleName()
176  << " off Element " << elm->GetName()
177  << " in " << mat->GetName()
178  << " Z= " << Z << " A= " << A
179  << " E(MeV)= " << part->GetKineticEnergy()/MeV << G4endl;
180  throw G4HadronicException(__FILE__, __LINE__,
181  " no applicable data set found for the isotope");
182  return 0.0;
183  //return dataSetList[idx]->ComputeCrossSection(part, elm, mat);
184 }
185 
186 G4double
187 G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* part,
188  G4int Z, G4int A,
189  const G4Isotope* iso,
190  const G4Element* elm,
191  const G4Material* mat)
192 {
193  for (G4int i = nDataSetList-1; i >= 0; --i) {
194  if (dataSetList[i]->IsIsoApplicable(part, Z, A, elm, mat) ) {
195  return dataSetList[i]->GetIsoCrossSection(part, Z, A, iso, elm, mat);
196  }
197  }
198  G4cout << "G4CrossSectionDataStore::GetCrossSection ERROR: "
199  << " no isotope cross section found"
200  << G4endl;
201  G4cout << " for " << part->GetDefinition()->GetParticleName()
202  << " off Element " << elm->GetName()
203  << " in " << mat->GetName()
204  << " Z= " << Z << " A= " << A
205  << " E(MeV)= " << part->GetKineticEnergy()/MeV << G4endl;
206  throw G4HadronicException(__FILE__, __LINE__,
207  " no applicable data set found for the isotope");
208  return 0.0;
209 }
210 
211 G4Element*
212 G4CrossSectionDataStore::SampleZandA(const G4DynamicParticle* part,
213  const G4Material* mat,
214  G4Nucleus& target)
215 {
216  G4int nElements = mat->GetNumberOfElements();
217  const G4ElementVector* theElementVector = mat->GetElementVector();
218  G4Element* anElement = (*theElementVector)[0];
219 
220  G4double cross = GetCrossSection(part, mat);
221 
222  // zero cross section case
223  if(0.0 >= cross) { return anElement; }
224 
225  // select element from a compound
226  if(1 < nElements) {
227  cross *= G4UniformRand();
228  for(G4int i=0; i<nElements; ++i) {
229  if(cross <= xsecelm[i]) {
230  anElement = (*theElementVector)[i];
231  break;
232  }
233  }
234  }
235 
236  G4int Z = G4lrint(anElement->GetZ());
237  G4Isotope* iso = 0;
238 
239  G4int i = nDataSetList-1;
240  if (dataSetList[i]->IsElementApplicable(part, Z, mat)) {
241 
242  //----------------------------------------------------------------
243  // element-wise cross section
244  // isotope cross section is not computed
245  //----------------------------------------------------------------
246  G4int nIso = anElement->GetNumberOfIsotopes();
247  if (0 >= nIso) {
248  G4cout << " Element " << anElement->GetName() << " Z= " << Z
249  << " has no isotopes " << G4endl;
250  throw G4HadronicException(__FILE__, __LINE__,
251  " Isotope vector is not defined");
252  return anElement;
253  }
254  // isotope abundances
255  G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
256  iso = (*isoVector)[0];
257 
258  // more than 1 isotope
259  if(1 < nIso) {
260  iso = dataSetList[i]->SelectIsotope(anElement, part->GetKineticEnergy());
261  }
262 
263  } else {
264 
265  //----------------------------------------------------------------
266  // isotope-wise cross section
267  // isotope cross section is computed
268  //----------------------------------------------------------------
269  G4int nIso = anElement->GetNumberOfIsotopes();
270  cross = 0.0;
271 
272  if (0 >= nIso) {
273  G4cout << " Element " << anElement->GetName() << " Z= " << Z
274  << " has no isotopes " << G4endl;
275  throw G4HadronicException(__FILE__, __LINE__,
276  " Isotope vector is not defined");
277  return anElement;
278  }
279 
280  // user-defined isotope abundances
281  G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
282  iso = (*isoVector)[0];
283 
284  // more than 1 isotope
285  if(1 < nIso) {
286  G4double* abundVector = anElement->GetRelativeAbundanceVector();
287  if(G4int(xseciso.size()) < nIso) { xseciso.resize(nIso); }
288 
289  for (G4int j = 0; j<nIso; ++j) {
290  G4double xsec = 0.0;
291  if(abundVector[j] > 0.0) {
292  iso = (*isoVector)[j];
293  xsec = abundVector[j]*
294  GetIsoCrossSection(part, Z, iso->GetN(), iso, anElement, mat, i);
295  }
296  cross += xsec;
297  xseciso[j] = cross;
298  }
299  cross *= G4UniformRand();
300  for (G4int j = 0; j<nIso; ++j) {
301  if(cross <= xseciso[j]) {
302  iso = (*isoVector)[j];
303  break;
304  }
305  }
306  }
307  }
308  target.SetIsotope(iso);
309  return anElement;
310 }
311 
312 void
313 G4CrossSectionDataStore::BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
314 {
315  if (nDataSetList == 0)
316  {
317  throw G4HadronicException(__FILE__, __LINE__,
318  "G4CrossSectionDataStore: no data sets registered");
319  return;
320  }
321  for (G4int i=0; i<nDataSetList; ++i) {
322  dataSetList[i]->BuildPhysicsTable(aParticleType);
323  }
324 }
325 
326 void
327 G4CrossSectionDataStore::DumpPhysicsTable(const G4ParticleDefinition& aParticleType)
328 {
329  // Print out all cross section data sets used and the energies at
330  // which they apply
331 
332  if (nDataSetList == 0) {
333  G4cout << "WARNING - G4CrossSectionDataStore::DumpPhysicsTable: "
334  << " no data sets registered" << G4endl;
335  return;
336  }
337 
338  for (G4int i = nDataSetList-1; i >= 0; --i) {
339  G4double e1 = dataSetList[i]->GetMinKinEnergy();
340  G4double e2 = dataSetList[i]->GetMaxKinEnergy();
341  if (i < nDataSetList-1) { G4cout << " "; }
342  G4cout << std::setw(25) << dataSetList[i]->GetName() << ": Emin(GeV)= "
343  << std::setw(4) << e1/GeV << " Emax(GeV)= "
344  << e2/GeV << G4endl;
345  if (dataSetList[i]->GetName() == "G4CrossSectionPairGG") {
346  dataSetList[i]->DumpPhysicsTable(aParticleType);
347  }
348  }
349 
350  G4cout << G4endl;
351 }
352 
353 void G4CrossSectionDataStore::DumpHtml(const G4ParticleDefinition&,
354  std::ofstream& outFile)
355 {
356  // Write cross section data set info to html physics list
357  // documentation page
358 
359  G4double ehi = 0;
360  G4double elo = 0;
361  for (G4int i = nDataSetList-1; i > 0; i--) {
362  elo = dataSetList[i]->GetMinKinEnergy()/GeV;
363  ehi = dataSetList[i]->GetMaxKinEnergy()/GeV;
364  outFile << " <li><b><a href=\"" << dataSetList[i]->GetName() << ".html\"> "
365  << dataSetList[i]->GetName() << "</a> from "
366  << elo << " GeV to " << ehi << " GeV </b></li>\n";
367  }
368 
369  G4double defaultHi = dataSetList[0]->GetMaxKinEnergy()/GeV;
370  if (ehi < defaultHi) {
371  outFile << " <li><b><a href=\"" << dataSetList[0]->GetName() << ".html\"> "
372  << dataSetList[0]->GetName() << "</a> from "
373  << ehi << " GeV to " << defaultHi << " GeV </b></li>\n";
374  }
375 }
Generated on Sat May 25 2013 14:33:41 for Geant4 by   doxygen 1.8.4