source: trunk/source/processes/electromagnetic/utils/src/G4EmSaturation.cc@ 1344

Last change on this file since 1344 was 1340, checked in by garnier, 15 years ago

update ti head
  • Property svn:executable set to *
File size: 8.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// $Id: G4EmSaturation.cc,v 1.11 2010/10/25 17:23:01 vnivanch Exp $
27// GEANT4 tag $Name: emutils-V09-03-23 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4EmSaturation
35//
36// Author: Vladimir Ivanchenko
37//
38// Creation date: 18.02.2008
39//
40// Modifications:
41//
42// -------------------------------------------------------------
43
44//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46
47#include "G4EmSaturation.hh"
48#include "G4LossTableManager.hh"
49#include "G4NistManager.hh"
50#include "G4Material.hh"
51#include "G4MaterialCutsCouple.hh"
52#include "G4Electron.hh"
53#include "G4Proton.hh"
54
55//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56
57G4EmSaturation::G4EmSaturation()
58{
59 verbose = 1;
60 manager = 0;
61
62 curMaterial = 0;
63 curBirks = 0.0;
64 curRatio = 1.0;
65 curChargeSq = 1.0;
66 nMaterials = 0;
67
68 electron = 0;
69 proton = 0;
70 nist = G4NistManager::Instance();
71
72 Initialise();
73}
74
75//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76
77G4EmSaturation::~G4EmSaturation()
78{}
79
80//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81
82G4double G4EmSaturation::VisibleEnergyDeposition(
83 const G4ParticleDefinition* p,
84 const G4MaterialCutsCouple* couple,
85 G4double length,
86 G4double edep,
87 G4double niel)
88{
89 if(edep <= 0.0) { return 0.0; }
90
91 G4double evis = edep;
92 G4double bfactor = FindBirksCoefficient(couple->GetMaterial());
93
94 if(bfactor > 0.0) {
95
96 G4int pdgCode = p->GetPDGEncoding();
97 // atomic relaxations for gamma incident
98 if(22 == pdgCode) {
99 evis /= (1.0 + bfactor*edep/manager->GetRange(electron,edep,couple));
100
101 // energy loss
102 } else {
103
104 // protections
105 G4double nloss = niel;
106 if(nloss < 0.0) nloss = 0.0;
107 G4double eloss = edep - nloss;
108
109 // neutrons
110 if(2112 == pdgCode || eloss < 0.0 || length <= 0.0) {
111 nloss = edep;
112 eloss = 0.0;
113 }
114
115 // continues energy loss
116 if(eloss > 0.0) { eloss /= (1.0 + bfactor*eloss/length); }
117
118 // non-ionizing energy loss
119 if(nloss > 0.0) {
120 if(!proton) { proton = G4Proton::Proton(); }
121 G4double escaled = nloss*curRatio;
122 G4double s = manager->GetRange(proton,escaled,couple)/curChargeSq;
123 nloss /= (1.0 + bfactor*nloss/s);
124 }
125
126 evis = eloss + nloss;
127 }
128 }
129
130 return evis;
131}
132
133//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
134
135G4double G4EmSaturation::FindG4BirksCoefficient(const G4Material* mat)
136{
137 G4String name = mat->GetName();
138 // is this material in the vector?
139
140 for(G4int j=0; j<nG4Birks; ++j) {
141 if(name == g4MatNames[j]) {
142 if(verbose > 0)
143 G4cout << "### G4EmSaturation::FindG4BirksCoefficient for "
144 << name << " is " << g4MatData[j]*MeV/mm << " mm/MeV "
145 << G4endl;
146 return g4MatData[j];
147 }
148 }
149 return FindBirksCoefficient(mat);
150}
151
152//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
153
154G4double G4EmSaturation::FindBirksCoefficient(const G4Material* mat)
155{
156 // electron should exist in any case
157 if(!manager) {
158 manager = G4LossTableManager::Instance();
159 electron= G4Electron::Electron();
160 }
161
162 if(mat == curMaterial) { return curBirks; }
163
164 curMaterial = mat;
165 curBirks = 0.0;
166 curRatio = 1.0;
167 curChargeSq = 1.0;
168
169 // seach in the run-time list
170 for(G4int i=0; i<nMaterials; ++i) {
171 if(mat == matPointers[i]) {
172 curBirks = mat->GetIonisation()->GetBirksConstant();
173 curRatio = massFactors[i];
174 curChargeSq = effCharges[i];
175 return curBirks;
176 }
177 }
178
179 G4String name = mat->GetName();
180 curBirks = mat->GetIonisation()->GetBirksConstant();
181
182 // material has no Birks coeffitient defined
183 // seach in the Geant4 list
184 if(curBirks == 0.0) {
185 for(G4int j=0; j<nG4Birks; ++j) {
186 if(name == g4MatNames[j]) {
187 mat->GetIonisation()->SetBirksConstant(g4MatData[j]);
188 curBirks = g4MatData[j];
189 break;
190 }
191 }
192 }
193
194 if(curBirks == 0.0 && verbose > 0) {
195 G4cout << "### G4EmSaturation::FindBirksCoefficient fails "
196 " for material " << name << G4endl;
197 }
198
199 // compute mean mass ratio
200 curRatio = 0.0;
201 curChargeSq = 0.0;
202 G4double norm = 0.0;
203 const G4ElementVector* theElementVector = mat->GetElementVector();
204 const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
205 size_t nelm = mat->GetNumberOfElements();
206 for (size_t i=0; i<nelm; ++i) {
207 const G4Element* elm = (*theElementVector)[i];
208 G4double Z = elm->GetZ();
209 G4double w = Z*Z*theAtomNumDensityVector[i];
210 curRatio += w/nist->GetAtomicMassAmu(G4int(Z));
211 curChargeSq = Z*Z*w;
212 norm += w;
213 }
214 curRatio *= proton_mass_c2/norm;
215 curChargeSq /= norm;
216
217 // store results
218 matPointers.push_back(mat);
219 matNames.push_back(name);
220 massFactors.push_back(curRatio);
221 effCharges.push_back(curChargeSq);
222 nMaterials++;
223 if(curBirks > 0.0 && verbose > 0) {
224 G4cout << "### G4EmSaturation::FindBirksCoefficient Birks coefficient for "
225 << name << " " << curBirks*MeV/mm << " mm/MeV" << G4endl;
226 }
227 return curBirks;
228}
229
230//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
231
232void G4EmSaturation::DumpBirksCoefficients()
233{
234 if(nMaterials > 0) {
235 G4cout << "### Birks coeffitients used in run time" << G4endl;
236 for(G4int i=0; i<nMaterials; ++i) {
237 G4double br = matPointers[i]->GetIonisation()->GetBirksConstant();
238 G4cout << " " << matNames[i] << " "
239 << br*MeV/mm << " mm/MeV" << " "
240 << br*matPointers[i]->GetDensity()*MeV*cm2/g
241 << " g/cm^2/MeV"
242 << G4endl;
243 }
244 }
245}
246
247//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
248
249void G4EmSaturation::DumpG4BirksCoefficients()
250{
251 if(nG4Birks > 0) {
252 G4cout << "### Birks coeffitients for Geant4 materials" << G4endl;
253 for(G4int i=0; i<nG4Birks; ++i) {
254 G4cout << " " << g4MatNames[i] << " "
255 << g4MatData[i]*MeV/mm << " mm/MeV" << G4endl;
256 }
257 }
258}
259
260//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
261
262void G4EmSaturation::Initialise()
263{
264 // M.Hirschberg et al., IEEE Trans. Nuc. Sci. 39 (1992) 511
265 // SCSN-38 kB = 0.00842 g/cm^2/MeV; rho = 1.06 g/cm^3
266 g4MatNames.push_back("G4_POLYSTYRENE");
267 g4MatData.push_back(0.07943*mm/MeV);
268
269 // C.Fabjan (private communication)
270 // kB = 0.006 g/cm^2/MeV; rho = 7.13 g/cm^3
271 g4MatNames.push_back("G4_BGO");
272 g4MatData.push_back(0.008415*mm/MeV);
273
274 // A.Ribon analysis of publications
275 // Scallettar et al., Phys. Rev. A25 (1982) 2419.
276 // NIM A 523 (2004) 275.
277 // kB = 0.022 g/cm^2/MeV; rho = 1.396 g/cm^3;
278 // ATLAS Efield = 10 kV/cm provide the strongest effect
279 g4MatNames.push_back("G4_lAr");
280 g4MatData.push_back(0.1576*mm/MeV);
281
282 //G4_BARIUM_FLUORIDE
283 //G4_CESIUM_IODIDE
284 //G4_GEL_PHOTO_EMULSION
285 //G4_PHOTO_EMULSION
286 //G4_PLASTIC_SC_VINYLTOLUENE
287 //G4_SODIUM_IODIDE
288 //G4_STILBENE
289 //G4_lAr
290 //G4_PbWO4
291 //G4_Lucite
292
293 nG4Birks = g4MatData.size();
294}
295
296//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Note: See TracBrowser for help on using the repository browser.