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