| Geant4 Cross Reference |
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 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4eIonisationCrossSectionHandler
33 //
34 // Author: V.Ivanchenko (Vladimir.Ivanchenko@cern.ch)
35 //
36 // Creation date: 25 Sept 2001
37 //
38 // Modifications:
39 // 10 Oct 2001 M.G. Pia Revision to improve code quality and consistency with design
40 // 19 Jul 2002 VI Create composite data set for material
41 // 21 Jan 2003 V.Ivanchenko Cut per region
42 // 28 Jan 2009 L.Pandola Added public method to make a easier migration of
43 // G4LowEnergyIonisation to G4LivermoreIonisationModel
44 // 15 Jul 2009 Nicolas A. Karakatsanis
45 //
46 // - BuildCrossSectionForMaterials method was revised in order to calculate the
47 // logarithmic values of the loaded data.
48 // It retrieves the data values from the G4EMLOW data files but, then, calculates the
49 // respective log values and loads them to seperate data structures.
50 // The EM data sets, initialized this way, contain both non-log and log values.
51 // These initialized data sets can enhance the computing performance of data interpolation
52 // operations
53 //
54 //
55 //
56 // -------------------------------------------------------------------
57
58 #include "G4eIonisationCrossSectionHandler.hh"
59 #include "G4SystemOfUnits.hh"
60 #include "G4VEnergySpectrum.hh"
61 #include "G4DataVector.hh"
62 #include "G4CompositeEMDataSet.hh"
63 #include "G4VDataSetAlgorithm.hh"
64 #include "G4LinLogLogInterpolation.hh"
65 #include "G4SemiLogInterpolation.hh"
66 #include "G4VEMDataSet.hh"
67 #include "G4EMDataSet.hh"
68 #include "G4Material.hh"
69 #include "G4ProductionCutsTable.hh"
70
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
72 G4eIonisationCrossSectionHandler::G4eIonisationCrossSectionHandler(
73 const G4VEnergySpectrum* spec, G4VDataSetAlgorithm* alg,
74 G4double emin, G4double emax, G4int nbin)
75 : G4VCrossSectionHandler(),
76 theParam(spec),verbose(0)
77 {
78 G4VCrossSectionHandler::Initialise(alg, emin, emax, nbin);
79 interp = new G4LinLogLogInterpolation();
80 }
81
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83
84 G4eIonisationCrossSectionHandler::~G4eIonisationCrossSectionHandler()
85 {
86 delete interp;
87 }
88
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90
91 std::vector<G4VEMDataSet*>* G4eIonisationCrossSectionHandler::BuildCrossSectionsForMaterials(
92 const G4DataVector& energyVector,
93 const G4DataVector* energyCuts)
94 {
95 std::vector<G4VEMDataSet*>* set = new std::vector<G4VEMDataSet*>;
96
97 G4DataVector* energies;
98 G4DataVector* cs;
99
100 G4DataVector* log_energies;
101 G4DataVector* log_cs;
102
103 std::size_t nOfBins = energyVector.size();
104
105 const G4ProductionCutsTable* theCoupleTable=
106 G4ProductionCutsTable::GetProductionCutsTable();
107 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
108
109 for (G4int mLocal=0; mLocal<numOfCouples; ++mLocal) {
110
111 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(mLocal);
112 const G4Material* material= couple->GetMaterial();
113 const G4ElementVector* elementVector = material->GetElementVector();
114 const G4double* nAtomsPerVolume = material->GetAtomicNumDensityVector();
115 G4int nElements = (G4int)material->GetNumberOfElements();
116
117 if(verbose > 0)
118 {
119 G4cout << "eIonisation CS for " << mLocal << "th material "
120 << material->GetName()
121 << " eEl= " << nElements << G4endl;
122 }
123
124 G4double tcut = (*energyCuts)[mLocal];
125
126 G4VDataSetAlgorithm* algo = interp->Clone();
127 G4VEMDataSet* setForMat = new G4CompositeEMDataSet(algo,1.,1.);
128
129 for (G4int i=0; i<nElements; ++i) {
130
131 G4int Z = (G4int) (*elementVector)[i]->GetZ();
132 G4int nShells = NumberOfComponents(Z);
133
134 energies = new G4DataVector;
135 cs = new G4DataVector;
136
137 log_energies = new G4DataVector;
138 log_cs = new G4DataVector;
139
140 G4double density = nAtomsPerVolume[i];
141
142 for (std::size_t bin=0; bin<nOfBins; ++bin) {
143
144 G4double e = energyVector[bin];
145 energies->push_back(e);
146 log_energies->push_back(std::log10(e));
147 G4double value = 0.0;
148 G4double log_value = -300;
149
150 if(e > tcut) {
151 for (G4int n=0; n<nShells; n++) {
152 G4double cross = FindValue(Z, e, n);
153 G4double p = theParam->Probability(Z, tcut, e, e, n);
154 value += cross * p * density;
155
156 if(verbose>0 && mLocal == 0 && e>=1. && e<=0.)
157 {
158 G4cout << "G4eIonCrossSH: e(MeV)= " << e/MeV
159 << " n= " << n
160 << " cross= " << cross
161 << " p= " << p
162 << " value= " << value
163 << " tcut(MeV)= " << tcut/MeV
164 << " rho= " << density
165 << " Z= " << Z
166 << G4endl;
167 }
168 }
169 if (value == 0.) value = 1e-300;
170 log_value = std::log10(value);
171 }
172 cs->push_back(value);
173 log_cs->push_back(log_value);
174 }
175 G4VDataSetAlgorithm* algoLocal = interp->Clone();
176 G4VEMDataSet* elSet = new G4EMDataSet(i,energies,cs,log_energies,log_cs,algoLocal,1.,1.);
177
178 setForMat->AddComponent(elSet);
179 }
180 set->push_back(setForMat);
181 }
182
183 return set;
184 }
185
186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
187
188 G4double G4eIonisationCrossSectionHandler::GetCrossSectionAboveThresholdForElement(G4double energy,
189 G4double cutEnergy,
190 G4int Z)
191 {
192 G4int nShells = NumberOfComponents(Z);
193 G4double value = 0.;
194 if(energy > cutEnergy)
195 {
196 for (G4int n=0; n<nShells; ++n) {
197 G4double cross = FindValue(Z, energy, n);
198 G4double p = theParam->Probability(Z, cutEnergy, energy, energy, n);
199 value += cross * p;
200 }
201 }
202 return value;
203 }
204