| Geant4 Cross Reference |
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 // GEANT4 Class file
29 //
30 // Description: Data on stopping power
31 //
32 // Author: Vladimir Ivanchenko
33 //
34 // Creation date: 23.10.2021
35 //
36 //----------------------------------------------------------------------------
37 //
38
39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
40
41 #include <fstream>
42 #include <sstream>
43
44 #include "G4IonICRU73Data.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4PhysicsLogVector.hh"
47 #include "G4PhysicsFreeVector.hh"
48 #include "G4EmParameters.hh"
49 #include "G4ProductionCutsTable.hh"
50 #include "G4MaterialCutsCouple.hh"
51 #include "G4Material.hh"
52 #include "G4Element.hh"
53
54 namespace {
55
56 const G4String namesICRU73[31] = {
57 "G4_A-150_TISSUE"
58 "G4_ADIPOSE_TISSUE_ICRP"
59 "G4_AIR"
60 "G4_ALUMINUM_OXIDE"
61 "G4_BONE_COMPACT_ICRU"
62 "G4_BONE_CORTICAL_ICRP"
63 "G4_C-552"
64 "G4_CALCIUM_FLUORIDE"
65 "G4_CARBON_DIOXIDE"
66 "G4_KAPTON"
67 "G4_LITHIUM_FLUORIDE"
68 "G4_LITHIUM_TETRABORATE"
69 "G4_LUCITE"
70 "G4_METHANE"
71 "G4_MUSCLE_STRIATED_ICRU"
72 "G4_MYLAR"
73 "G4_NYLON-6-6"
74 "G4_PHOTO_EMULSION"
75 "G4_PLASTIC_SC_VINYLTOLUENE"
76 "G4_POLYCARBONATE"
77 "G4_POLYETHYLENE"
78 "G4_POLYSTYRENE"
79 "G4_PROPANE"
80 "G4_Pyrex_Glass"
81 "G4_SILICON_DIOXIDE"
82 "G4_SODIUM_IODIDE"
83 "G4_TEFLON"
84 "G4_TISSUE-METHANE"
85 "G4_TISSUE-PROPANE"
86 "G4_WATER"
87 "G4_WATER_VAPOR"};
88 const G4String namesICRU90[3] = {"G4_AIR","G4_GRAPHITE","G4_WATER"};
89 const G4double densityCoef[3] = {0.996, 1.025, 0.998};
90 const G4int NZ = 27;
91 const G4int zdat[NZ] = { 5, 6, 7, 8, 13, 14, 15, 16, 22, 26,
92 28, 29, 30, 31, 32, 33, 34, 40, 47, 48,
93 49, 51, 52, 72, 73, 74, 79 };
94 }
95
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
97
98 G4IonICRU73Data::G4IonICRU73Data()
99 {
100 fEmin = 0.025*CLHEP::MeV;
101 fEmax = 2.5*CLHEP::MeV;
102 fNbins = fNbinsPerDecade*G4lrint(std::log10(fEmax/fEmin));
103 fVector = new G4PhysicsFreeVector(fSpline);
104 for(G4int i=3; i<=ZPROJMAX; ++i) {
105 fMatData[i] = new std::vector<G4PhysicsLogVector*>;
106 }
107 }
108
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
110
111 G4IonICRU73Data::~G4IonICRU73Data()
112 {
113 delete fVector;
114 for(G4int i=3; i<=ZPROJMAX; ++i) {
115 auto v = fMatData[i];
116 if(nullptr != v) {
117 for(auto & dat : *v) {
118 delete dat;
119 }
120 delete v;
121 }
122 for(G4int j=1; j<=ZTARGMAX; ++j) {
123 delete fElmData[i][j];
124 }
125 }
126 }
127
128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
129
130 G4double G4IonICRU73Data::GetDEDX(const G4Material* mat, const G4int Z,
131 const G4double e, const G4double loge) const
132 {
133 // if data does not exist the result is zero
134 if (Z > ZPROJMAX) { return 0.0; }
135
136 G4PhysicsLogVector* v = nullptr;
137 if (1 == mat->GetNumberOfElements()) {
138 G4int Z1 = (*(mat->GetElementVector()))[0]->GetZasInt();
139 if(Z1 > ZTARGMAX) { return 0.0; }
140 v = fElmData[Z][Z1];
141 }
142 else {
143 G4int idx = fMatIndex[mat->GetIndex()];
144 if (idx < 0) { return 0.0; }
145 v = (*(fMatData[Z]))[idx];
146 }
147 if(nullptr == v) { return 0.0; }
148 G4double res = (e > fEmin) ? v->LogVectorValue(e, loge)
149 : (*v)[0]*std::sqrt(e/fEmin);
150 return res;
151 }
152
153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
154
155 void G4IonICRU73Data::Initialise()
156 {
157 // fill directory path
158 if(fDataDirectory.empty()) {
159 std::ostringstream ost;
160 ost << G4EmParameters::Instance()->GetDirLEDATA() << "/ion_stopping_data/";
161 fDataDirectory = ost.str();
162 }
163
164 const std::size_t nmat = G4Material::GetNumberOfMaterials();
165
166 // do nothing if for a new run list of materials is not changed
167 if(nmat == fMatIndex.size()) { return; }
168
169 // look over all materials and recompute all materials
170 // do not recompute element data if exist
171 if(1 < fVerbose) {
172 G4cout << "### G4IonICRU73Data::Initialise() for " << nmat
173 << " materials" << G4endl;
174 }
175 fMatIndex.resize(nmat, -1);
176 for(G4int j=3; j<=ZPROJMAX; ++j) {
177 fMatData[j]->resize(nmat, nullptr);
178 }
179 G4bool useICRU90 = G4EmParameters::Instance()->UseICRU90Data();
180 auto mtable = G4Material::GetMaterialTable();
181
182 for(G4int i=0; i<(G4int)nmat; ++i) {
183 const G4Material* mat = (*mtable)[i];
184 const G4String matname = mat->GetName();
185 G4int idx = (G4int)mat->GetIndex();
186 if(1 < fVerbose) {
187 G4cout << i << ". material: " << matname
188 << " idx=" << idx << " matIdx=" << fMatIndex[idx]
189 << G4endl;
190 }
191 if(fMatIndex[idx] == -1) {
192 G4bool isOK = false;
193
194 // for material in ICRU90
195 if(useICRU90) {
196 for(G4int j=0; j<3; ++j) {
197 if(matname == namesICRU90[j]) {
198 ReadMaterialData(mat, densityCoef[j], true);
199 isOK = true;
200 if(1 < fVerbose) {
201 G4cout << "ICRU90 material " << matname << G4endl;
202 }
203 break;
204 }
205 }
206 }
207 // for material in ICRU73
208 if(!isOK) {
209 for(G4int j=0; j<31; ++j) {
210 if(matname == namesICRU73[j]) {
211 ReadMaterialData(mat, 1.0, false);
212 isOK = true;
213 if(1 < fVerbose) {
214 G4cout << "ICRU73 material " << matname << G4endl;
215 }
216 break;
217 }
218 }
219 }
220 // dEdx is combined from element stopping power
221 // Target element Z <= ZTARGMAX
222 if(!isOK) {
223 const std::size_t nelm = mat->GetNumberOfElements();
224 auto elmv = mat->GetElementVector();
225 G4bool isOut = false;
226 for (std::size_t j=0; j<nelm; ++j) {
227 if ((*elmv)[j]->GetZasInt() > ZTARGMAX) {
228 isOut = true;
229 break;
230 }
231 }
232 if (!isOut) {
233 ReadElementData(mat, useICRU90);
234 isOK = true;
235 if(1 < fVerbose) {
236 G4cout << "Data via elements for " << matname << G4endl;
237 }
238 }
239 }
240 if(isOK) { fMatIndex[idx] = i; }
241 }
242 if(1 < fVerbose) {
243 G4cout << " matData: " << fMatData[i] << G4endl;
244 }
245 }
246 }
247
248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
249
250 void G4IonICRU73Data::ReadMaterialData(const G4Material* mat,
251 const G4double coeff,
252 const G4bool useICRU90)
253 {
254 const G4String name = mat->GetName();
255 const G4int idx = (G4int)mat->GetIndex();
256 for(G4int Z=3; Z<=ZPROJMAX; ++Z) {
257 std::ostringstream ost;
258 ost << fDataDirectory << "icru";
259 G4int Z1 = Z;
260 G4double scale = 1.0;
261 if(useICRU90 && Z <= 18) {
262 ost << "90";
263 } else {
264 ost << "73";
265 for(G4int i=0; i<NZ; ++i) {
266 if(Z == zdat[i]) {
267 break;
268 } else if(i == NZ-1) {
269 Z1 = zdat[NZ - 1];
270 scale = (G4double)(Z*Z)/(G4double)(Z1*Z1);
271 } else if(Z > zdat[i] && Z < zdat[i+1]) {
272 if(Z - zdat[i] <= zdat[i + 1] - Z) {
273 Z1 = zdat[i];
274 } else {
275 Z1 = zdat[i+1];
276 }
277 scale = (G4double)(Z*Z)/(G4double)(Z1*Z1);
278 break;
279 }
280 }
281 }
282 if(nullptr == (*(fMatData[Z1]))[idx]) {
283 ost << "/z" << Z1 << "_" << name << ".dat";
284 G4PhysicsLogVector* v = RetrieveVector(ost, false);
285 if(nullptr != v) {
286 const G4double fact = coeff *
287 mat->GetDensity() * CLHEP::MeV * 1000 * CLHEP::cm2 / CLHEP::g;
288 v->ScaleVector(CLHEP::MeV, fact);
289 if(2 < fVerbose) {
290 G4cout << "### Data for " << name
291 << " and projectile Z=" << Z1 << G4endl;
292 G4cout << *v << G4endl;
293 }
294 (*(fMatData[Z1]))[idx] = v;
295 }
296 }
297 if(Z != Z1) {
298 auto v2 = (*(fMatData[Z1]))[idx];
299 if(nullptr != v2) {
300 auto v1 = new G4PhysicsLogVector(*v2);
301 (*(fMatData[Z]))[idx] = v1;
302 v1->ScaleVector(1.0, scale);
303 }
304 }
305 }
306 }
307
308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
309
310 void G4IonICRU73Data::ReadElementData(const G4Material* mat, G4bool useICRU90)
311 {
312 const G4ElementVector* elmv = mat->GetElementVector();
313 const G4double* dens = mat->GetFractionVector();
314 const G4int nelm = (G4int)mat->GetNumberOfElements();
315 for(G4int Z=3; Z<=ZPROJMAX; ++Z) {
316 if (1 < fVerbose) {
317 G4cout << "ReadElementData for " << mat->GetName() << " Z=" << Z
318 << " Nelm=" << nelm << G4endl;
319 }
320 G4PhysicsLogVector* v = nullptr;
321 if(1 == nelm) {
322 v = FindOrBuildElementData(Z, (*elmv)[0]->GetZasInt(), useICRU90);
323 } else {
324 v = new G4PhysicsLogVector(fEmin, fEmax, fNbins, fSpline);
325 for(G4int i=0; i<=fNbins; ++i) {
326 G4double dedx = 0.0;
327 for(G4int j=0; j<nelm; ++j) {
328 G4PhysicsLogVector* v1 =
329 FindOrBuildElementData(Z, (*elmv)[j]->GetZasInt(), useICRU90);
330 dedx += (*v1)[i]*dens[j];
331 }
332 v->PutValue(i, dedx);
333 }
334 if(fSpline) { v->FillSecondDerivatives(); }
335 (*(fMatData[Z]))[(G4int)mat->GetIndex()] = v;
336 }
337 // scale data for correct units
338 if(nullptr != v) {
339 const G4double fact =
340 mat->GetDensity() * CLHEP::MeV * 1000 * CLHEP::cm2 / CLHEP::g;
341 v->ScaleVector(CLHEP::MeV, fact);
342 if(2 < fVerbose) {
343 G4cout << "### Data for "<< mat->GetName()
344 << " for projectile Z=" << Z << G4endl;
345 G4cout << *v << G4endl;
346 }
347 }
348 }
349 }
350
351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
352
353 G4PhysicsLogVector*
354 G4IonICRU73Data::FindOrBuildElementData(const G4int Z, const G4int Z1,
355 G4bool useICRU90)
356 {
357 G4PhysicsLogVector* v = nullptr;
358 if(Z <= ZPROJMAX && Z1 <= ZTARGMAX) {
359 v = fElmData[Z][Z1];
360 if(nullptr == v) {
361 G4int Z2 = Z1;
362 G4bool isICRU90 = (useICRU90 && Z <= 18 &&
363 (Z1 == 1 || Z1 == 6 || Z1 == 7 || Z1 == 8));
364
365 G4double scale = 1.0;
366 if(!isICRU90) {
367 for(G4int i=0; i<NZ; ++i) {
368 if(Z1 == zdat[i]) {
369 break;
370 } else if(i == NZ-1) {
371 Z2 = zdat[NZ - 1];
372 scale = (G4double)Z1/(G4double)Z2;
373 } else if(Z1 > zdat[i] && Z1 < zdat[i+1]) {
374 if(Z1 - zdat[i] <= zdat[i + 1] - Z1) {
375 Z2 = zdat[i];
376 } else {
377 Z2 = zdat[i+1];
378 }
379 scale = (G4double)Z1/(G4double)Z2;
380 break;
381 }
382 }
383 }
384
385 std::ostringstream ost;
386 ost << fDataDirectory << "icru";
387 if(isICRU90) { ost << "90"; }
388 else { ost << "73"; }
389 ost << "/z" << Z << "_" << Z2 << ".dat";
390 v = RetrieveVector(ost, false);
391 fElmData[Z][Z2] = v;
392 if(Z1 != Z2 && nullptr != v) {
393 fElmData[Z][Z1] = new G4PhysicsLogVector(*v);
394 fElmData[Z][Z1]->ScaleVector(1.0, scale);
395 }
396 }
397 }
398 return v;
399 }
400
401 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
402
403 G4PhysicsLogVector*
404 G4IonICRU73Data::RetrieveVector(std::ostringstream& ost, G4bool warn)
405 {
406 G4PhysicsLogVector* v = nullptr;
407 std::ifstream filein(ost.str().c_str());
408 if (!filein.is_open()) {
409 if(warn) {
410 G4ExceptionDescription ed;
411 ed << "Data file <" << ost.str().c_str()
412 << "> is not opened";
413 G4Exception("G4IonICRU73Data::RetrieveVector(..)","em013",
414 FatalException, ed, "Check G4LEDATA");
415 }
416 } else {
417 if(fVerbose > 0) {
418 G4cout << "File " << ost.str()
419 << " is opened by G4IonICRU73Data" << G4endl;
420 }
421
422 // retrieve data from DB
423 if(!fVector->Retrieve(filein, true)) {
424 G4ExceptionDescription ed;
425 ed << "Data file <" << ost.str().c_str()
426 << "> is not retrieved!";
427 G4Exception("G4IonICRU73Data::RetrieveVector(..)","had015",
428 FatalException, ed, "Check G4LEDATA");
429 } else {
430 if(fSpline) { fVector->FillSecondDerivatives(); }
431 fVector->EnableLogBinSearch(G4EmParameters::Instance()->NumberForFreeVector());
432 v = new G4PhysicsLogVector(fEmin, fEmax, fNbins, fSpline);
433 for(G4int i=0; i<=fNbins; ++i) {
434 G4double e = v->Energy(i);
435 G4double dedx = fVector->Value(e);
436 v->PutValue(i, dedx);
437 }
438 if(fSpline) { v->FillSecondDerivatives(); }
439 if(fVerbose > 2) { G4cout << *v << G4endl; }
440 }
441 }
442 return v;
443 }
444
445 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
446