Geant4  10.01.p03
G4EmSaturation.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: G4EmSaturation.cc 81936 2014-06-06 15:42:55Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4EmSaturation
34 //
35 // Author: Vladimir Ivanchenko
36 //
37 // Creation date: 18.02.2008
38 //
39 // Modifications:
40 //
41 // -------------------------------------------------------------
42 
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 
46 #include "G4EmSaturation.hh"
47 #include "G4PhysicalConstants.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "G4LossTableManager.hh"
50 #include "G4NistManager.hh"
51 #include "G4Material.hh"
52 #include "G4MaterialCutsCouple.hh"
53 #include "G4ParticleTable.hh"
54 
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56 
58  : manager(0)
59 {
60  verbose = verb;
61  manager = 0;
62 
63  curMaterial = 0;
64  curBirks = 0.0;
65  curRatio = 1.0;
66  curChargeSq = 1.0;
67  nMaterials = nWarnings = 0;
68 
69  electron = 0;
70  proton = 0;
72 
74 }
75 
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77 
79 {}
80 
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
82 
84  const G4ParticleDefinition* p,
85  const G4MaterialCutsCouple* couple,
86  G4double length,
87  G4double edep,
88  G4double niel)
89 {
90  if(edep <= 0.0) { return 0.0; }
91 
92  G4double evis = edep;
93  G4double bfactor = FindBirksCoefficient(couple->GetMaterial());
94 
95  if(bfactor > 0.0) {
96 
97  G4int pdgCode = p->GetPDGEncoding();
98  // atomic relaxations for gamma incident
99  if(22 == pdgCode && electron) {
100  //G4cout << "%% gamma edep= " << edep/keV << " keV " <<manager << G4endl;
101  evis /= (1.0 + bfactor*edep/manager->GetRange(electron,edep,couple));
102 
103  // energy loss
104  } else {
105 
106  // protections
107  G4double nloss = niel;
108  if(nloss < 0.0) { nloss = 0.0; }
109  G4double eloss = edep - nloss;
110 
111  // neutrons and neutral hadrons
112  if(0.0 == p->GetPDGCharge() || eloss < 0.0 || length <= 0.0) {
113  nloss = edep;
114  eloss = 0.0;
115  }
116 
117  // continues energy loss
118  if(eloss > 0.0) { eloss /= (1.0 + bfactor*eloss/length); }
119 
120  // non-ionizing energy loss
121  if(nloss > 0.0 && proton) {
122  G4double escaled = nloss*curRatio;
123  /*
124  G4cout << "%% p edep= " << nloss/keV << " keV Escaled= "
125  << escaled << " MeV in " << couple->GetMaterial()->GetName()
126  << " " << p->GetParticleName()
127  << G4endl;
128  */
129  G4double range = manager->GetRange(proton,escaled,couple)/curChargeSq;
130  nloss /= (1.0 + bfactor*nloss/range);
131  }
132 
133  evis = eloss + nloss;
134  }
135  }
136 
137  return evis;
138 }
139 
140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
141 
143 {
144  G4String name = mat->GetName();
145  // is this material in the vector?
146 
147  for(G4int j=0; j<nG4Birks; ++j) {
148  if(name == g4MatNames[j]) {
149  if(verbose > 0)
150  G4cout << "### G4EmSaturation::FindG4BirksCoefficient for "
151  << name << " is " << g4MatData[j]*MeV/mm << " mm/MeV "
152  << G4endl;
153  return g4MatData[j];
154  }
155  }
156  return FindBirksCoefficient(mat);
157 }
158 
159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
160 
162 {
163  // electron and proton should exist in any case
164  if(!manager) {
168  }
169 
170  curMaterial = mat;
171  curBirks = 0.0;
172  curRatio = 1.0;
173  curChargeSq = 1.0;
174 
175  // seach in the run-time list
176  for(G4int i=0; i<nMaterials; ++i) {
177  if(mat == matPointers[i]) {
179  curRatio = massFactors[i];
180  curChargeSq = effCharges[i];
181  return;
182  }
183  }
184 
185  G4String name = mat->GetName();
187 
188  // material has no Birks coeffitient defined
189  // seach in the Geant4 list
190  if(curBirks == 0.0) {
191  for(G4int j=0; j<nG4Birks; ++j) {
192  if(name == g4MatNames[j]) {
194  curBirks = g4MatData[j];
195  break;
196  }
197  }
198  }
199 
200  if(curBirks == 0.0) {
201  if(0 < nWarnings) {
202  ++nWarnings;
204  ed << "Birks constants are not defined for material " << name
205  << " ! \n Define Birks constants for the material"
206  << " or not apply saturation.";
207  G4Exception("G4EmSaturation::InitialiseBirksCoefficient", "em0088",
208  JustWarning, ed);
209  }
210  return;
211  }
212 
213  // compute mean mass ratio
214  curRatio = 0.0;
215  curChargeSq = 0.0;
216  G4double norm = 0.0;
217  const G4ElementVector* theElementVector = mat->GetElementVector();
218  const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
219  size_t nelm = mat->GetNumberOfElements();
220  for (size_t i=0; i<nelm; ++i) {
221  const G4Element* elm = (*theElementVector)[i];
222  G4double Z = elm->GetZ();
223  G4double w = Z*Z*theAtomNumDensityVector[i];
225  curChargeSq = Z*Z*w;
226  norm += w;
227  }
228  curRatio *= proton_mass_c2/norm;
229  curChargeSq /= norm;
230 
231  // store results
232  matPointers.push_back(mat);
233  matNames.push_back(name);
234  massFactors.push_back(curRatio);
235  effCharges.push_back(curChargeSq);
236  nMaterials++;
237  if(verbose > 0) {
238  G4cout << "### G4EmSaturation::FindBirksCoefficient Birks coefficient for "
239  << name << " " << curBirks*MeV/mm << " mm/MeV" << G4endl;
240  }
241  return;
242 }
243 
244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
245 
247 {
248  if(nMaterials > 0) {
249  G4cout << "### Birks coeffitients used in run time" << G4endl;
250  for(G4int i=0; i<nMaterials; ++i) {
251  G4double br = matPointers[i]->GetIonisation()->GetBirksConstant();
252  G4cout << " " << matNames[i] << " "
253  << br*MeV/mm << " mm/MeV" << " "
254  << br*matPointers[i]->GetDensity()*MeV*cm2/g
255  << " g/cm^2/MeV"
256  << G4endl;
257  }
258  }
259 }
260 
261 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
262 
264 {
265  if(nG4Birks > 0) {
266  G4cout << "### Birks coeffitients for Geant4 materials" << G4endl;
267  for(G4int i=0; i<nG4Birks; ++i) {
268  G4cout << " " << g4MatNames[i] << " "
269  << g4MatData[i]*MeV/mm << " mm/MeV" << G4endl;
270  }
271  }
272 }
273 
274 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
275 
277 {
278  // M.Hirschberg et al., IEEE Trans. Nuc. Sci. 39 (1992) 511
279  // SCSN-38 kB = 0.00842 g/cm^2/MeV; rho = 1.06 g/cm^3
280  g4MatNames.push_back("G4_POLYSTYRENE");
281  g4MatData.push_back(0.07943*mm/MeV);
282 
283  // C.Fabjan (private communication)
284  // kB = 0.006 g/cm^2/MeV; rho = 7.13 g/cm^3
285  g4MatNames.push_back("G4_BGO");
286  g4MatData.push_back(0.008415*mm/MeV);
287 
288  // A.Ribon analysis of publications
289  // Scallettar et al., Phys. Rev. A25 (1982) 2419.
290  // NIM A 523 (2004) 275.
291  // kB = 0.022 g/cm^2/MeV; rho = 1.396 g/cm^3;
292  // ATLAS Efield = 10 kV/cm provide the strongest effect
293  g4MatNames.push_back("G4_lAr");
294  g4MatData.push_back(0.1576*mm/MeV);
295 
296  //G4_BARIUM_FLUORIDE
297  //G4_CESIUM_IODIDE
298  //G4_GEL_PHOTO_EMULSION
299  //G4_PHOTO_EMULSION
300  //G4_PLASTIC_SC_VINYLTOLUENE
301  //G4_SODIUM_IODIDE
302  //G4_STILBENE
303  //G4_lAr
304  //G4_PbWO4
305  //G4_Lucite
306 
307  nG4Birks = g4MatData.size();
308 }
309 
310 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
std::vector< const G4Material * > matPointers
const G4ParticleDefinition * electron
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:226
static const double MeV
Definition: G4SIunits.hh:193
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
static G4LossTableManager * Instance()
static const double cm2
Definition: G4SIunits.hh:107
std::vector< G4double > g4MatData
const G4ParticleDefinition * proton
std::vector< G4Element * > G4ElementVector
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double FindBirksCoefficient(const G4Material *)
void InitialiseG4materials()
G4NistManager * nist
G4String name
Definition: TRTMaterials.hh:40
const G4String & GetName() const
Definition: G4Material.hh:178
void InitialiseBirksCoefficient(const G4Material *)
virtual G4double VisibleEnergyDeposition(const G4ParticleDefinition *, const G4MaterialCutsCouple *, G4double length, G4double edepTotal, G4double edepNIEL=0.0)
void SetBirksConstant(G4double value)
virtual ~G4EmSaturation()
const G4Material * curMaterial
G4double GetZ() const
Definition: G4Element.hh:131
std::vector< G4double > effCharges
std::vector< G4String > g4MatNames
void DumpBirksCoefficients()
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
G4double GetBirksConstant() const
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:206
std::vector< G4String > matNames
G4GLOB_DLL std::ostream G4cout
G4double FindG4BirksCoefficient(const G4Material *)
std::vector< G4double > massFactors
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void DumpG4BirksCoefficients()
static G4ParticleTable * GetParticleTable()
static const double g
Definition: G4SIunits.hh:162
G4LossTableManager * manager
G4double GetAtomicMassAmu(const G4String &symb) const
#define G4endl
Definition: G4ios.hh:61
G4double curChargeSq
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
G4EmSaturation(G4int verb)
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
G4double GetRange(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
static const double mm
Definition: G4SIunits.hh:102
const G4Material * GetMaterial() const