Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 105120 2017-07-13 13:24:43Z 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 
57 G4int G4EmSaturation::nMaterials = 0;
58 std::vector<G4double> G4EmSaturation::massFactors;
59 std::vector<G4double> G4EmSaturation::effCharges;
60 std::vector<G4double> G4EmSaturation::g4MatData;
61 std::vector<G4String> G4EmSaturation::g4MatNames;
62 
64 {
65  verbose = verb;
66 
67  nWarnings = nG4Birks = 0;
68 
69  electron = nullptr;
70  proton = nullptr;
71  nist = G4NistManager::Instance();
72 }
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 
77 {}
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80 
82  const G4ParticleDefinition* p,
83  const G4MaterialCutsCouple* couple,
84  G4double length,
85  G4double edep,
86  G4double niel) const
87 {
88  if(edep <= 0.0) { return 0.0; }
89 
90  G4double evis = edep;
91  G4double bfactor = couple->GetMaterial()->GetIonisation()->GetBirksConstant();
92 
93  if(bfactor > 0.0) {
94 
95  // atomic relaxations for gamma incident
96  if(22 == p->GetPDGEncoding()) {
97  //G4cout << "%% gamma edep= " << edep/keV << " keV " <<manager << G4endl;
98  evis /= (1.0 + bfactor*edep/
99  G4LossTableManager::Instance()->GetRange(electron,edep,couple));
100 
101  // energy loss
102  } else {
103 
104  // protections
105  G4double nloss = std::max(niel, 0.0);
106  G4double eloss = edep - nloss;
107 
108  // neutrons and neutral hadrons
109  if(0.0 == p->GetPDGCharge() || eloss < 0.0 || length <= 0.0) {
110  nloss = edep;
111  eloss = 0.0;
112  } else {
113 
114  // continues energy loss
115  eloss /= (1.0 + bfactor*eloss/length);
116  }
117  // non-ionizing energy loss
118  if(nloss > 0.0) {
119  G4int idx = couple->GetMaterial()->GetIndex();
120  G4double escaled = nloss*massFactors[idx];
121  /*
122  G4cout << "%% p edep= " << nloss/keV << " keV Escaled= "
123  << escaled << " MeV in " << couple->GetMaterial()->GetName()
124  << " " << p->GetParticleName()
125  << G4endl;
126  */
128  ->GetRange(proton,escaled,couple)/effCharges[idx];
129  nloss /= (1.0 + bfactor*nloss/range);
130  }
131  evis = eloss + nloss;
132  }
133  }
134  return evis;
135 }
136 
137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
138 
140 {
141  nMaterials = G4Material::GetNumberOfMaterials();
142  massFactors.resize(nMaterials, 1.0);
143  effCharges.resize(nMaterials, 1.0);
144 
145  if(0 == nG4Birks) { InitialiseG4materials(); }
146 
147  for(G4int i=0; i<nMaterials; ++i) {
148  InitialiseBirksCoefficient((*G4Material::GetMaterialTable())[i]);
149  }
150  if(verbose > 0) { DumpBirksCoefficients(); }
151 }
152 
153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
154 
156 {
157  if(0 == nG4Birks) { InitialiseG4materials(); }
158 
159  G4String name = mat->GetName();
160  // is this material in the vector?
161 
162  for(G4int j=0; j<nG4Birks; ++j) {
163  if(name == g4MatNames[j]) {
164  if(verbose > 0)
165  G4cout << "### G4EmSaturation::FindG4BirksCoefficient for "
166  << name << " is " << g4MatData[j]*MeV/mm << " mm/MeV "
167  << G4endl;
168  return g4MatData[j];
169  }
170  }
171  return 0.0;
172 }
173 
174 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
175 
176 void G4EmSaturation::InitialiseBirksCoefficient(const G4Material* mat)
177 {
178  // electron and proton should exist in any case
179  if(!electron) {
181  proton = G4ParticleTable::GetParticleTable()->FindParticle("proton");
182  if(!electron || !proton) {
183  G4Exception("G4EmSaturation::InitialiseBirksCoefficient", "em0001",
184  FatalException, "both electron and proton should exist");
185  }
186  }
187 
188  G4double curBirks = mat->GetIonisation()->GetBirksConstant();
189 
190  G4String name = mat->GetName();
191 
192  // material has no Birks coeffitient defined
193  // seach in the Geant4 list
194  if(curBirks == 0.0) {
195  for(G4int j=0; j<nG4Birks; ++j) {
196  if(name == g4MatNames[j]) {
197  mat->GetIonisation()->SetBirksConstant(g4MatData[j]);
198  curBirks = g4MatData[j];
199  break;
200  }
201  }
202  }
203 
204  if(curBirks == 0.0) { return; }
205 
206  // compute mean mass ratio
207  G4double curRatio = 0.0;
208  G4double curChargeSq = 0.0;
209  G4double norm = 0.0;
210  const G4ElementVector* theElementVector = mat->GetElementVector();
211  const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
212  size_t nelm = mat->GetNumberOfElements();
213  for (size_t i=0; i<nelm; ++i) {
214  const G4Element* elm = (*theElementVector)[i];
215  G4double Z = elm->GetZ();
216  G4double w = Z*Z*theAtomNumDensityVector[i];
217  curRatio += w/nist->GetAtomicMassAmu(G4int(Z));
218  curChargeSq = Z*Z*w;
219  norm += w;
220  }
221  curRatio *= proton_mass_c2/norm;
222  curChargeSq /= norm;
223 
224  // store results
225  G4int idx = mat->GetIndex();
226  massFactors[idx] = curRatio;
227  effCharges[idx] = curChargeSq;
228 }
229 
230 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
231 
233 {
234  G4cout << "### Birks coefficients used in run time" << G4endl;
236  for(G4int i=0; i<nMaterials; ++i) {
237  const G4Material* mat = (*mtable)[i];
238  G4double br = mat->GetIonisation()->GetBirksConstant();
239  if(br > 0.0) {
240  G4cout << " " << mat->GetName() << " "
241  << br*MeV/mm << " mm/MeV" << " "
242  << br*mat->GetDensity()*MeV*cm2/g
243  << " g/cm^2/MeV massFactor= " << massFactors[i]
244  << " effCharge= " << effCharges[i] << G4endl;
245  }
246  }
247 }
248 
249 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
250 
252 {
253  if(nG4Birks > 0) {
254  G4cout << "### Birks coefficients for Geant4 materials" << G4endl;
255  for(G4int i=0; i<nG4Birks; ++i) {
256  G4cout << " " << g4MatNames[i] << " "
257  << g4MatData[i]*MeV/mm << " mm/MeV" << G4endl;
258  }
259  }
260 }
261 
262 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
263 
264 void G4EmSaturation::InitialiseG4materials()
265 {
266  nG4Birks = 4;
267  g4MatData.reserve(nG4Birks);
268 
269  // M.Hirschberg et al., IEEE Trans. Nuc. Sci. 39 (1992) 511
270  // SCSN-38 kB = 0.00842 g/cm^2/MeV; rho = 1.06 g/cm^3
271  g4MatNames.push_back("G4_POLYSTYRENE");
272  g4MatData.push_back(0.07943*mm/MeV);
273 
274  // C.Fabjan (private communication)
275  // kB = 0.006 g/cm^2/MeV; rho = 7.13 g/cm^3
276  g4MatNames.push_back("G4_BGO");
277  g4MatData.push_back(0.008415*mm/MeV);
278 
279  // A.Ribon analysis of publications
280  // Scallettar et al., Phys. Rev. A25 (1982) 2419.
281  // NIM A 523 (2004) 275.
282  // kB = 0.022 g/cm^2/MeV; rho = 1.396 g/cm^3;
283  // ATLAS Efield = 10 kV/cm provide the strongest effect
284  // kB = 0.1576*mm/MeV
285  // A. Kiryunin and P.Strizenec "Geant4 hadronic
286  // working group meeting " kB = 0.041/9.13 g/cm^2/MeV
287  g4MatNames.push_back("G4_lAr");
288  g4MatData.push_back(0.032*mm/MeV);
289 
290  //G4_BARIUM_FLUORIDE
291  //G4_CESIUM_IODIDE
292  //G4_GEL_PHOTO_EMULSION
293  //G4_PHOTO_EMULSION
294  //G4_PLASTIC_SC_VINYLTOLUENE
295  //G4_SODIUM_IODIDE
296  //G4_STILBENE
297  //G4_lAr
298 
299  //G4_PbWO4 - CMS value
300  g4MatNames.push_back("G4_PbWO4");
301  g4MatData.push_back(0.0333333*mm/MeV);
302 
303  //G4_Lucite
304 
305 }
306 
307 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
const XML_Char * name
Definition: expat.h:151
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:226
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
static G4LossTableManager * Instance()
static constexpr double mm
Definition: G4SIunits.hh:115
std::vector< G4Element * > G4ElementVector
static constexpr double cm2
Definition: G4SIunits.hh:120
size_t GetIndex() const
Definition: G4Material.hh:262
const char * p
Definition: xmltok.h:285
const G4String & GetName() const
Definition: G4Material.hh:178
void SetBirksConstant(G4double value)
virtual ~G4EmSaturation()
G4double GetZ() const
Definition: G4Element.hh:131
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:587
std::vector< G4Material * > G4MaterialTable
G4double GetDensity() const
Definition: G4Material.hh:180
virtual G4double VisibleEnergyDeposition(const G4ParticleDefinition *, const G4MaterialCutsCouple *, G4double length, G4double edepTotal, G4double edepNIEL=0.0) const
void DumpBirksCoefficients()
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
G4double GetBirksConstant() const
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:206
G4GLOB_DLL std::ostream G4cout
const G4ParticleDefinition const G4Material *G4double range
void InitialiseG4Saturation()
G4double FindG4BirksCoefficient(const G4Material *)
float proton_mass_c2
Definition: hepunit.py:275
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:594
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void DumpG4BirksCoefficients()
static G4ParticleTable * GetParticleTable()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double GetAtomicMassAmu(const G4String &symb) const
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
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)
const G4Material * GetMaterial() const