Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PWATotalXsecTable.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 implementation file
31 //
32 // File name: G4PWATotalXsecTable
33 //
34 // Author: Mihaly Novak
35 //
36 // Creation date: 18.05.2015
37 //
38 // Class description:
39 // Class to load and handle elastic, first and second transport cross sections
40 // precomputed by using ELSEPA [1] in the 100 eV - 1 GeV kinetic and Z = 1-103
41 // energy range for electrons and positrons.G4PWATotalXsecZ is responsible to
42 // to handle cross sections by individual Z that are used in the current
43 // geometry and G4PWATotalXsecTable is a collection of G4PWATotalXsecZ objects.
44 //
45 // Modifications:
46 //
47 // References:
48 // [1] Francesc Salvat, Aleksander Jablonski, Cedric J Powell,
49 // ELSEPA—Dirac partial-wave calculation of elastic scattering of electrons
50 // and positrons by atoms, positive ions and molecules,
51 // Computer physics communications; 165, 2, (2005)
52 //
53 // -----------------------------------------------------------------------------
54 
55 
56 #include "G4PWATotalXsecTable.hh"
57 
58 
59 #include <vector>
60 #include <iostream>
61 #include <fstream>
62 #include <cstdlib>
63 #include <cmath>
64 
65 #include "G4MaterialTable.hh"
66 #include "G4Material.hh"
67 #include "G4Log.hh"
68 #include "G4Exp.hh"
69 
71 // G4PWATotalXsecZ: sub-class for PWA xsec data that belong to a given Z number
73 
75 const G4double G4PWATotalXsecZ::fgPWATotalXsecEnergyGrid[]={
76  // energy bin values for total elastic, first and second transport cross scetions in MeV
77  1.00000000e-04, 1.16591440e-04, 1.35935639e-04, 1.58489319e-04, 1.84784980e-04, 2.15443469e-04, 2.51188643e-04, 2.92864456e-04,
78  3.41454887e-04, 3.98107171e-04, 4.64158883e-04, 5.41169527e-04, 6.30957344e-04, 7.35642254e-04, 8.57695899e-04, 1.00000000e-03,
79  1.16591440e-03, 1.35935639e-03, 1.58489319e-03, 1.84784980e-03, 2.15443469e-03, 2.51188643e-03, 2.92864456e-03, 3.41454887e-03,
80  3.98107171e-03, 4.64158883e-03, 5.41169527e-03, 6.30957344e-03, 7.35642254e-03, 8.57695899e-03, 1.00000000e-02, 1.16591440e-02,
81  1.35935639e-02, 1.58489319e-02, 1.84784980e-02, 2.15443469e-02, 2.51188643e-02, 2.92864456e-02, 3.41454887e-02, 3.98107171e-02,
82  4.64158883e-02, 5.41169527e-02, 6.30957344e-02, 7.35642254e-02, 8.57695899e-02, 1.00000000e-01, 1.16591440e-01, 1.35935639e-01,
83  1.58489319e-01, 1.84784980e-01, 2.15443469e-01, 2.51188643e-01, 2.92864456e-01, 3.41454887e-01, 3.98107171e-01, 4.64158883e-01,
84  5.41169527e-01, 6.30957344e-01, 7.35642254e-01, 8.57695899e-01, 1.00000000e+00, 1.16591440e+00, 1.35935639e+00, 1.58489319e+00,
85  1.84784980e+00, 2.15443469e+00, 2.51188643e+00, 2.92864456e+00, 3.41454887e+00, 3.98107171e+00, 4.64158883e+00, 5.41169527e+00,
86  6.30957344e+00, 7.35642254e+00, 8.57695899e+00, 1.00000000e+01, 1.16591440e+01, 1.35935639e+01, 1.58489319e+01, 1.84784980e+01,
87  2.15443469e+01, 2.51188643e+01, 2.92864456e+01, 3.41454887e+01, 3.98107171e+01, 4.64158883e+01, 5.41169527e+01, 6.30957344e+01,
88  7.35642254e+01, 8.57695899e+01, 1.00000000e+02, 1.16591440e+02, 1.35935639e+02, 1.58489319e+02, 1.84784980e+02, 2.15443469e+02,
89  2.51188643e+02, 2.92864456e+02, 3.41454887e+02, 3.98107171e+02, 4.64158883e+02, 5.41169527e+02, 6.30957344e+02, 7.35642254e+02,
90  8.57695899e+02, 1.00000000e+03
91  };
92 
94 G4PWATotalXsecZ::G4PWATotalXsecZ(G4int Z){
95  G4int nn = fgNumTotalXsecBins*6;
96  for(G4int i=0; i<nn; ++i) {
97  fPWAXsecs[i] = 0.0;
98  fInterpParamA[i] = 0.0;
99  fInterpParamB[i] = 0.0;
100  }
101  LoadPWATotalXsecZ(Z);
102 }
103 
105 void G4PWATotalXsecZ::LoadPWATotalXsecZ(G4int Z){
106  G4double dum;
107  char fname[512];
108  char* path = getenv("G4LEDATA");
109  if (!path) {
110  G4Exception("G4PWATotalXsecZ::LoadPWATotalXsecZ()","em0006",
112  "Environment variable G4LEDATA not defined");
113  return;
114  }
115 
116  std::string pathString(path);
117  sprintf(fname,"%s/msc_GS/xsecs/xsecs_%d",path,Z);
118  std::ifstream infile(fname,std::ios::in);
119  if(!infile.is_open()){
120  char msgc[512];
121  sprintf(msgc," Total PWA xsection %s not found.",fname);
122  G4Exception("G4PWATotalXsecZ::LoadPWATotalXsecZ()","em0006",
124  msgc);
125  return;
126  }
127  G4double dummy;
128 
129  for(G4int i=0; i<fgNumTotalXsecBins; ++i)
130  for(G4int j=0; j<7; ++j)
131  if(j==0) infile >> dum;
132  else {
133  // load pwa xsection that are stored in cm2 units in file and change to
134  // Geant4 internal length2 units
135  infile >> dummy;
136  fPWAXsecs[(j-1)*fgNumTotalXsecBins+i] = dummy*CLHEP::cm2;
137  }
138  infile.close();
139 
140  // compute log-log linear intrp. parameters
141  for(G4int i=0; i<fgNumTotalXsecBins-1; ++i)
142  for(G4int k=0; k<6; ++k) {
143  G4int j = k*fgNumTotalXsecBins+i;
144  G4double val2 = fPWAXsecs[j+1];
145  G4double val1 = fPWAXsecs[j];
146 
147  fInterpParamA[j] = G4Log(val2/val1)/G4Log(fgPWATotalXsecEnergyGrid[i+1]/fgPWATotalXsecEnergyGrid[i]);
148  fInterpParamB[j] = G4Exp(G4Log(val1) - fInterpParamA[j]*G4Log(fgPWATotalXsecEnergyGrid[i]));
149  }
150 }
151 
153 // Get the index of the lower energy bin edge
155  // log(fgPWATotalXsecEnergyGrid[0]);
156  const G4double lne0 = -9.21034037197618e+00;
157  // 1./log(fgPWATotalXsecEnergyGrid[i+1]/fgPWATotalXsecEnergyGrid[i]);
158  const G4double invlnde = 6.51441722854880e+00;
159 
160  return (G4int)((G4Log(energy)-lne0)*invlnde);
161 }
162 
164 // j-dependent type interploated cross section in Geant4 internal length2 unit
165 // energy is assumed to be in [MeV]
167  // protection : out of energy grid range
168  if(energy < GetLowestEnergy())
169  return GetLowestXsecValue(j);
170  if(energy >= GetHighestEnergy())
171  return GetHighestXsecValue(j);
172 
173  // normal case log-log linear intrp.
174  G4int k = j*fgNumTotalXsecBins+elowindex;
175  return G4Exp(G4Log(energy)*fInterpParamA[k])*fInterpParamB[k];
176 }
177 
180  // protection : out of energy grid range
181  if(energy < GetLowestEnergy())
182  return GetLowestXsecValue(j);
183  if(energy >= GetHighestEnergy())
184  return GetHighestXsecValue(j);
185 
186  // normal case log-log linear intrp.
187  G4int elowindex = GetPWATotalXsecEnergyBinIndex(energy);
188  G4int k = j*fgNumTotalXsecBins+elowindex;
189  return G4Exp(G4Log(energy)*fInterpParamA[k])*fInterpParamB[k];
190 }
191 
192 
194 // G4PWATotalXsecTable
196 
197 G4PWATotalXsecZ* G4PWATotalXsecTable::fgPWATotalXsecTable[fgNumZet] = {0};
198 
201  for(G4int i = 0; i < fgNumZet; ++i)
202  if(fgPWATotalXsecTable[i]) {
203  delete fgPWATotalXsecTable[i];
204  fgPWATotalXsecTable[i] = 0;
205  }
206 }
207 
210  G4int isUsedZ[fgNumZet] ={0}; //xsec data available up to fgNumZet Z-number
211  // check used elements
212  G4MaterialTable *theMaterialTable = G4Material::GetMaterialTable();
213  for(unsigned int imat = 0; imat < theMaterialTable->size(); ++imat) {
214  const G4ElementVector *theElemVect = ((*theMaterialTable)[imat])->GetElementVector();
215  for(unsigned int ielem = 0; ielem < theElemVect->size(); ++ielem) {
216  G4int zet = G4lrint((*theElemVect)[ielem]->GetZ());
217  zet = zet>fgNumZet ? fgNumZet : zet;
218  if(!isUsedZ[zet-1])
219  isUsedZ[zet-1] = 1;
220  }
221  }
222 
223  for(G4int i = 0; i < fgNumZet; ++i)
224  if(isUsedZ[i] && !fgPWATotalXsecTable[i]) // used but not there yet -> load it
225  fgPWATotalXsecTable[i] = new G4PWATotalXsecZ(i+1);
226  else if(!isUsedZ[i] && fgPWATotalXsecTable[i]) { // there but not used now -> delete
227  delete fgPWATotalXsecTable[i];
228  fgPWATotalXsecTable[i] = 0;
229  }
230 }
231 
std::vector< G4Element * > G4ElementVector
G4double GetHighestEnergy() const
G4double GetLowestEnergy() const
static constexpr double cm2
G4double GetInterpXsec(G4double energy, G4int elowindex, G4int j) const
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:587
std::vector< G4Material * > G4MaterialTable
int G4int
Definition: G4Types.hh:78
G4int GetPWATotalXsecEnergyBinIndex(G4double energy) const
G4double GetLowestXsecValue(G4int j) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
int G4lrint(double ad)
Definition: templates.hh:163
G4double energy(const ThreeVector &p, const G4double m)
string fname
Definition: test.py:308
G4double GetHighestXsecValue(G4int j) const
double G4double
Definition: G4Types.hh:76