Geant4  10.02.p01
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 92921 2015-09-21 15:06:51Z 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(nullptr)
59 {
60  verbose = verb;
61 
62  curMaterial = nullptr;
63  curBirks = 0.0;
64  curRatio = 1.0;
65  curChargeSq = 1.0;
66  nMaterials = nWarnings = 0;
67 
68  electron = nullptr;
69  proton = nullptr;
71 
73 }
74 
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76 
78 {}
79 
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81 
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 && electron) {
99  //G4cout << "%% gamma edep= " << edep/keV << " keV " <<manager << G4endl;
100  evis /= (1.0 + bfactor*edep/manager->GetRange(electron,edep,couple));
101 
102  // energy loss
103  } else {
104 
105  // protections
106  G4double nloss = niel;
107  if(nloss < 0.0) { nloss = 0.0; }
108  G4double eloss = edep - nloss;
109 
110  // neutrons and neutral hadrons
111  if(0.0 == p->GetPDGCharge() || eloss < 0.0 || length <= 0.0) {
112  nloss = edep;
113  eloss = 0.0;
114  }
115 
116  // continues energy loss
117  if(eloss > 0.0) { eloss /= (1.0 + bfactor*eloss/length); }
118 
119  // non-ionizing energy loss
120  if(nloss > 0.0 && proton) {
121  G4double escaled = nloss*curRatio;
122  /*
123  G4cout << "%% p edep= " << nloss/keV << " keV Escaled= "
124  << escaled << " MeV in " << couple->GetMaterial()->GetName()
125  << " " << p->GetParticleName()
126  << G4endl;
127  */
128  G4double range = manager->GetRange(proton,escaled,couple)/curChargeSq;
129  nloss /= (1.0 + bfactor*nloss/range);
130  }
131 
132  evis = eloss + nloss;
133  }
134  }
135 
136  return evis;
137 }
138 
139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
140 
142 {
143  G4String name = mat->GetName();
144  // is this material in the vector?
145 
146  for(G4int j=0; j<nG4Birks; ++j) {
147  if(name == g4MatNames[j]) {
148  if(verbose > 0)
149  G4cout << "### G4EmSaturation::FindG4BirksCoefficient for "
150  << name << " is " << g4MatData[j]*MeV/mm << " mm/MeV "
151  << G4endl;
152  return g4MatData[j];
153  }
154  }
155  return FindBirksCoefficient(mat);
156 }
157 
158 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
159 
161 {
162  // electron and proton should exist in any case
163  if(!manager) {
167  }
168 
169  curMaterial = mat;
170  curBirks = 0.0;
171  curRatio = 1.0;
172  curChargeSq = 1.0;
173 
174  // seach in the run-time list
175  for(G4int i=0; i<nMaterials; ++i) {
176  if(mat == matPointers[i]) {
178  curRatio = massFactors[i];
179  curChargeSq = effCharges[i];
180  return;
181  }
182  }
183 
184  G4String name = mat->GetName();
186 
187  // material has no Birks coeffitient defined
188  // seach in the Geant4 list
189  if(curBirks == 0.0) {
190  for(G4int j=0; j<nG4Birks; ++j) {
191  if(name == g4MatNames[j]) {
193  curBirks = g4MatData[j];
194  break;
195  }
196  }
197  }
198 
199  if(curBirks == 0.0) {
200  if(0 < nWarnings) {
201  ++nWarnings;
203  ed << "Birks constants are not defined for material " << name
204  << " ! \n Define Birks constants for the material"
205  << " or not apply saturation.";
206  G4Exception("G4EmSaturation::InitialiseBirksCoefficient", "em0088",
207  JustWarning, ed);
208  }
209  return;
210  }
211 
212  // compute mean mass ratio
213  curRatio = 0.0;
214  curChargeSq = 0.0;
215  G4double norm = 0.0;
216  const G4ElementVector* theElementVector = mat->GetElementVector();
217  const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
218  size_t nelm = mat->GetNumberOfElements();
219  for (size_t i=0; i<nelm; ++i) {
220  const G4Element* elm = (*theElementVector)[i];
221  G4double Z = elm->GetZ();
222  G4double w = Z*Z*theAtomNumDensityVector[i];
224  curChargeSq = Z*Z*w;
225  norm += w;
226  }
227  curRatio *= proton_mass_c2/norm;
228  curChargeSq /= norm;
229 
230  // store results
231  matPointers.push_back(mat);
232  matNames.push_back(name);
233  massFactors.push_back(curRatio);
234  effCharges.push_back(curChargeSq);
235  nMaterials++;
236  if(verbose > 0) {
237  G4cout << "### G4EmSaturation::FindBirksCoefficient Birks coefficient for "
238  << name << " " << curBirks*MeV/mm << " mm/MeV" << G4endl;
239  }
240  return;
241 }
242 
243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
244 
246 {
247  if(nMaterials > 0) {
248  G4cout << "### Birks coeffitients used in run time" << G4endl;
249  for(G4int i=0; i<nMaterials; ++i) {
250  G4double br = matPointers[i]->GetIonisation()->GetBirksConstant();
251  G4cout << " " << matNames[i] << " "
252  << br*MeV/mm << " mm/MeV" << " "
253  << br*matPointers[i]->GetDensity()*MeV*cm2/g
254  << " g/cm^2/MeV"
255  << G4endl;
256  }
257  }
258 }
259 
260 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
261 
263 {
264  if(nG4Birks > 0) {
265  G4cout << "### Birks coeffitients for Geant4 materials" << G4endl;
266  for(G4int i=0; i<nG4Birks; ++i) {
267  G4cout << " " << g4MatNames[i] << " "
268  << g4MatData[i]*MeV/mm << " mm/MeV" << G4endl;
269  }
270  }
271 }
272 
273 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
274 
276 {
277  // M.Hirschberg et al., IEEE Trans. Nuc. Sci. 39 (1992) 511
278  // SCSN-38 kB = 0.00842 g/cm^2/MeV; rho = 1.06 g/cm^3
279  g4MatNames.push_back("G4_POLYSTYRENE");
280  g4MatData.push_back(0.07943*mm/MeV);
281 
282  // C.Fabjan (private communication)
283  // kB = 0.006 g/cm^2/MeV; rho = 7.13 g/cm^3
284  g4MatNames.push_back("G4_BGO");
285  g4MatData.push_back(0.008415*mm/MeV);
286 
287  // A.Ribon analysis of publications
288  // Scallettar et al., Phys. Rev. A25 (1982) 2419.
289  // NIM A 523 (2004) 275.
290  // kB = 0.022 g/cm^2/MeV; rho = 1.396 g/cm^3;
291  // ATLAS Efield = 10 kV/cm provide the strongest effect
292  g4MatNames.push_back("G4_lAr");
293  g4MatData.push_back(0.1576*mm/MeV);
294 
295  //G4_BARIUM_FLUORIDE
296  //G4_CESIUM_IODIDE
297  //G4_GEL_PHOTO_EMULSION
298  //G4_PHOTO_EMULSION
299  //G4_PLASTIC_SC_VINYLTOLUENE
300  //G4_SODIUM_IODIDE
301  //G4_STILBENE
302  //G4_lAr
303  //G4_PbWO4
304  //G4_Lucite
305 
306  nG4Birks = g4MatData.size();
307 }
308 
309 //....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:211
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
static G4LossTableManager * Instance()
static const double cm2
Definition: G4SIunits.hh:119
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
const G4double w[NPOINTSGL]
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:180
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:114
const G4Material * GetMaterial() const