| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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........oo
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","G
89 const G4double densityCoef[3] = {0.996, 1.02
90 const G4int NZ = 27;
91 const G4int zdat[NZ] = { 5, 6, 7, 8, 13,
92 28, 29, 30, 31, 32,
93 49, 51, 52, 72, 73,
94 }
95
96 //....oooOO0OOooo........oooOO0OOooo........oo
97
98 G4IonICRU73Data::G4IonICRU73Data()
99 {
100 fEmin = 0.025*CLHEP::MeV;
101 fEmax = 2.5*CLHEP::MeV;
102 fNbins = fNbinsPerDecade*G4lrint(std::log10(
103 fVector = new G4PhysicsFreeVector(fSpline);
104 for(G4int i=3; i<=ZPROJMAX; ++i) {
105 fMatData[i] = new std::vector<G4PhysicsLog
106 }
107 }
108
109 //....oooOO0OOooo........oooOO0OOooo........oo
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........oo
129
130 G4double G4IonICRU73Data::GetDEDX(const G4Mate
131 const G4doub
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]
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->LogVectorVal
149 : (*v)[0]*std::sqrt(e/fEmin);
150 return res;
151 }
152
153 //....oooOO0OOooo........oooOO0OOooo........oo
154
155 void G4IonICRU73Data::Initialise()
156 {
157 // fill directory path
158 if(fDataDirectory.empty()) {
159 std::ostringstream ost;
160 ost << G4EmParameters::Instance()->GetDirL
161 fDataDirectory = ost.str();
162 }
163
164 const std::size_t nmat = G4Material::GetNumb
165
166 // do nothing if for a new run list of mater
167 if(nmat == fMatIndex.size()) { return; }
168
169 // look over all materials and recompute all
170 // do not recompute element data if exist
171 if(1 < fVerbose) {
172 G4cout << "### G4IonICRU73Data::Initialise
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(
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: " << matnam
188 << " idx=" << idx << " matIdx="
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[
199 isOK = true;
200 if(1 < fVerbose) {
201 G4cout << "ICRU90 material " <<
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 " <<
215 }
216 break;
217 }
218 }
219 }
220 // dEdx is combined from element stoppin
221 // Target element Z <= ZTARGMAX
222 if(!isOK) {
223 const std::size_t nelm = mat->GetNumberOfEle
224 auto elmv = mat->GetElementVector();
225 G4bool isOut = false;
226 for (std::size_t j=0; j<nelm; ++j) {
227 if ((*elmv)[j]->GetZasInt() > ZTARGM
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 " << ma
237 }
238 }
239 }
240 if(isOK) { fMatIndex[idx] = i; }
241 }
242 if(1 < fVerbose) {
243 G4cout << " matData: " << fMatData[i
244 }
245 }
246 }
247
248 //....oooOO0OOooo........oooOO0OOooo........oo
249
250 void G4IonICRU73Data::ReadMaterialData(const G
251 const G
252 const G
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 << ".da
284 G4PhysicsLogVector* v = RetrieveVector(o
285 if(nullptr != v) {
286 const G4double fact = coeff *
287 mat->GetDensity() * CLHEP::MeV * 1000 * CL
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........oo
309
310 void G4IonICRU73Data::ReadElementData(const G4
311 {
312 const G4ElementVector* elmv = mat->GetElemen
313 const G4double* dens = mat->GetFractionVecto
314 const G4int nelm = (G4int)mat->GetNumberOfEl
315 for(G4int Z=3; Z<=ZPROJMAX; ++Z) {
316 if (1 < fVerbose) {
317 G4cout << "ReadElementData for " << mat-
318 << " Nelm=" << nelm << G4endl;
319 }
320 G4PhysicsLogVector* v = nullptr;
321 if(1 == nelm) {
322 v = FindOrBuildElementData(Z, (*elmv)[0]
323 } else {
324 v = new G4PhysicsLogVector(fEmin, fEmax,
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)[
330 dedx += (*v1)[i]*dens[j];
331 }
332 v->PutValue(i, dedx);
333 }
334 if(fSpline) { v->FillSecondDerivatives()
335 (*(fMatData[Z]))[(G4int)mat->GetIndex()]
336 }
337 // scale data for correct units
338 if(nullptr != v) {
339 const G4double fact =
340 mat->GetDensity() * CLHEP::MeV * 1000
341 v->ScaleVector(CLHEP::MeV, fact);
342 if(2 < fVerbose) {
343 G4cout << "### Data for "<< mat->GetNa
344 << " for projectile Z=" << Z <<
345 G4cout << *v << G4endl;
346 }
347 }
348 }
349 }
350
351 //....oooOO0OOooo........oooOO0OOooo........oo
352
353 G4PhysicsLogVector*
354 G4IonICRU73Data::FindOrBuildElementData(const
355 G4bool
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 |
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] - Z
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........oo
402
403 G4PhysicsLogVector*
404 G4IonICRU73Data::RetrieveVector(std::ostringst
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::RetrieveVe
414 FatalException, ed, "Check G
415 }
416 } else {
417 if(fVerbose > 0) {
418 G4cout << "File " << ost.str()
419 << " is opened by G4IonICRU73Data
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::RetrieveVe
428 FatalException, ed, "Check G
429 } else {
430 if(fSpline) { fVector->FillSecondDerivat
431 fVector->EnableLogBinSearch(G4EmParamete
432 v = new G4PhysicsLogVector(fEmin, fEmax,
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 << G4end
440 }
441 }
442 return v;
443 }
444
445 //....oooOO0OOooo........oooOO0OOooo........oo
446