| 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 // 10.06.97 created. V. Grichine
27 // 18.11.98 simplified public interface; new methods for materials. mma
28 // 31.01.01 redesign of ComputeMatSandiaMatrix(). mma
29 // 16.02.01 adapted for STL. mma
30 // 22.02.01 GetsandiaCofForMaterial(energy) return 0 below lowest interval mma
31 // 03.04.01 fnulcof returned if energy < emin
32 // 10.07.01 Migration to STL. M. Verderi.
33 // 03.02.04 Update distructor V.Ivanchenko
34 // 05.03.04 New methods for old sorting algorithm for PAI model. V.Grichine
35 // 26.10.11 new scheme for G4Exception (mma)
36 // 22.05.13 preparation of material table without dynamical arrays. V. Grichine
37 // 09.07.14 modify low limit in GetSandiaCofPerAtom() and material. mma
38 // 10.07.14 modify low limit for water. VI
39
40 #include "G4SandiaTable.hh"
41
42 #include "G4Material.hh"
43 #include "G4MaterialTable.hh"
44 #include "G4PhysicalConstants.hh"
45 #include "G4StaticSandiaData.hh"
46 #include "G4SystemOfUnits.hh"
47
48 const G4double G4SandiaTable::funitc[5] = {CLHEP::keV, CLHEP::cm2* CLHEP::keV / CLHEP::g,
49 CLHEP::cm2* CLHEP::keV* CLHEP::keV / CLHEP::g,
50 CLHEP::cm2* CLHEP::keV* CLHEP::keV* CLHEP::keV / CLHEP::g,
51 CLHEP::cm2* CLHEP::keV* CLHEP::keV* CLHEP::keV* CLHEP::keV / CLHEP::g};
52
53 G4int G4SandiaTable::fCumulInterval[] = {0};
54
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
56
57 G4SandiaTable::G4SandiaTable(const G4Material* material) : fMaterial(material)
58 {
59 fMatSandiaMatrix = nullptr;
60 fMatSandiaMatrixPAI = nullptr;
61 fPhotoAbsorptionCof = nullptr;
62
63 fMatNbOfIntervals = 0;
64
65 fMaxInterval = 0;
66 fVerbose = 0;
67
68 // build the CumulInterval array
69 if (0 == fCumulInterval[0]) {
70 fCumulInterval[0] = 1;
71
72 for (G4int Z = 1; Z < 101; ++Z) {
73 fCumulInterval[Z] = fCumulInterval[Z - 1] + fNbOfIntervals[Z];
74 }
75 }
76
77 fMaxInterval = 0;
78 fSandiaCofPerAtom.resize(4, 0.0);
79 fLowerI1 = false;
80 // compute macroscopic Sandia coefs for a material
81 ComputeMatSandiaMatrix(); // mma
82 }
83
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
85
86 G4SandiaTable::~G4SandiaTable()
87 {
88 if (fMatSandiaMatrix != nullptr) {
89 fMatSandiaMatrix->clearAndDestroy();
90 delete fMatSandiaMatrix;
91 }
92 if (fMatSandiaMatrixPAI != nullptr) {
93 fMatSandiaMatrixPAI->clearAndDestroy();
94 delete fMatSandiaMatrixPAI;
95 }
96
97 delete[] fPhotoAbsorptionCof;
98 }
99
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
101
102 void G4SandiaTable::GetSandiaCofPerAtom(
103 G4int Z, G4double energy, std::vector<G4double>& coeff) const
104 {
105 #ifdef G4VERBOSE
106 if (Z < 1 || Z > 100) {
107 Z = PrintErrorZ(Z, "GetSandiaCofPerAtom");
108 }
109 if (4 > coeff.size()) {
110 PrintErrorV("GetSandiaCofPerAtom(): input vector is resized");
111 coeff.resize(4);
112 }
113 #endif
114 G4double Emin = fSandiaTable[fCumulInterval[Z - 1]][0] * CLHEP::keV;
115 // G4double Iopot = fIonizationPotentials[Z]*eV;
116 // if (Emin < Iopot) Emin = Iopot;
117
118 G4int row = 0;
119 if (energy <= Emin) {
120 energy = Emin;
121 }
122 else {
123 G4int interval = fNbOfIntervals[Z] - 1;
124 row = fCumulInterval[Z - 1] + interval;
125 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
126 while ((interval > 0) && (energy < fSandiaTable[row][0] * CLHEP::keV)) {
127 --interval;
128 row = fCumulInterval[Z - 1] + interval;
129 }
130 }
131
132 G4double AoverAvo = Z * amu / fZtoAratio[Z];
133
134 coeff[0] = AoverAvo * funitc[1] * fSandiaTable[row][1];
135 coeff[1] = AoverAvo * funitc[2] * fSandiaTable[row][2];
136 coeff[2] = AoverAvo * funitc[3] * fSandiaTable[row][3];
137 coeff[3] = AoverAvo * funitc[4] * fSandiaTable[row][4];
138 }
139
140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
141
142 void G4SandiaTable::GetSandiaCofWater(G4double energy, std::vector<G4double>& coeff) const
143 {
144 #ifdef G4VERBOSE
145 if (4 > coeff.size()) {
146 PrintErrorV("GetSandiaCofWater: input vector is resized");
147 coeff.resize(4);
148 }
149 #endif
150 G4int i = 0;
151 if (energy > fH2OlowerI1[0][0] * CLHEP::keV) {
152 i = fH2OlowerInt - 1;
153 for (; i > 0; --i) {
154 if (energy >= fH2OlowerI1[i][0] * CLHEP::keV) {
155 break;
156 }
157 }
158 }
159 coeff[0] = funitc[1] * fH2OlowerI1[i][1];
160 coeff[1] = funitc[2] * fH2OlowerI1[i][2];
161 coeff[2] = funitc[3] * fH2OlowerI1[i][3];
162 coeff[3] = funitc[4] * fH2OlowerI1[i][4];
163 }
164
165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
166
167 G4double G4SandiaTable::GetWaterEnergyLimit() const
168 {
169 return fH2OlowerI1[fH2OlowerInt - 1][0] * CLHEP::keV;
170 }
171
172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
173
174 G4double G4SandiaTable::GetWaterCofForMaterial(G4int i, G4int j) const
175 {
176 return fH2OlowerI1[i][j] * funitc[j];
177 }
178
179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
180
181 G4double G4SandiaTable::GetZtoA(G4int Z)
182 {
183 #ifdef G4VERBOSE
184 if (Z < 1 || Z > 100) {
185 Z = PrintErrorZ(Z, "GetSandiaCofPerAtom");
186 }
187 #endif
188 return fZtoAratio[Z];
189 }
190
191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
192
193 #ifdef G4VERBOSE
194
195 G4int G4SandiaTable::PrintErrorZ(G4int Z, const G4String& ss)
196 {
197 G4String sss = "G4SandiaTable::" + ss + "()";
198 G4ExceptionDescription ed;
199 ed << "Atomic number out of range Z= " << Z << "; closest value is used";
200 G4Exception(sss, "mat060", JustWarning, ed, "");
201 return (Z > 100) ? 100 : 1;
202 }
203
204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
205
206 void G4SandiaTable::PrintErrorV(const G4String& ss)
207 {
208 G4String sss = "G4SandiaTable::" + ss;
209 G4ExceptionDescription ed;
210 G4Exception(sss, "mat061", JustWarning, "Wrong input parameters");
211 }
212 #endif
213
214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
215
216 void G4SandiaTable::ComputeMatSandiaMatrix()
217 {
218 // get list of elements
219 const auto NbElm = (G4int)fMaterial->GetNumberOfElements();
220 const G4ElementVector* ElementVector = fMaterial->GetElementVector();
221
222 auto Z = new G4int[NbElm]; // Atomic number
223
224 // determine the maximum number of energy-intervals for this material
225 G4int MaxIntervals = 0;
226 G4int elm, z;
227
228 // here we compute only for a mixture, so no waring or exception
229 // if z is out of validity interval
230 for (elm = 0; elm < NbElm; ++elm) {
231 z = G4lrint((*ElementVector)[elm]->GetZ());
232 if (z < 1) {
233 z = 1;
234 }
235 else if (z > 100) {
236 z = 100;
237 }
238 Z[elm] = z;
239 MaxIntervals += fNbOfIntervals[z];
240 }
241
242 // copy the Energy bins in a tmp1 array
243 //(take care of the Ionization Potential of each element)
244 auto tmp1 = new G4double[MaxIntervals];
245 G4double IonizationPot;
246 G4int interval1 = 0;
247
248 for (elm = 0; elm < NbElm; ++elm) {
249 z = Z[elm];
250 IonizationPot = fIonizationPotentials[z] * CLHEP::eV;
251 for (G4int row = fCumulInterval[z - 1]; row < fCumulInterval[z]; ++row) {
252 tmp1[interval1] = std::max(fSandiaTable[row][0] * CLHEP::keV, IonizationPot);
253 ++interval1;
254 }
255 }
256 // sort the energies in strickly increasing values in a tmp2 array
257 //(eliminate redondances)
258
259 auto tmp2 = new G4double[MaxIntervals];
260 G4double Emin;
261 G4int interval2 = 0;
262
263 do {
264 Emin = DBL_MAX;
265
266 for (G4int i1 = 0; i1 < MaxIntervals; ++i1) {
267 Emin = std::min(Emin, tmp1[i1]); // find the minimum
268 }
269 if (Emin < DBL_MAX) {
270 tmp2[interval2] = Emin;
271 ++interval2;
272 }
273 // copy Emin in tmp2
274 for (G4int j1 = 0; j1 < MaxIntervals; ++j1) {
275 if (tmp1[j1] <= Emin) {
276 tmp1[j1] = DBL_MAX;
277 } // eliminate from tmp1
278 }
279 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
280 } while (Emin < DBL_MAX);
281
282 // create the sandia matrix for this material
283
284 fMatSandiaMatrix = new G4OrderedTable();
285 G4int interval;
286
287 for (interval = 0; interval < interval2; ++interval) {
288 fMatSandiaMatrix->push_back(new G4DataVector(5, 0.));
289 }
290
291 // ready to compute the Sandia coefs for the material
292
293 const G4double* NbOfAtomsPerVolume = fMaterial->GetVecNbOfAtomsPerVolume();
294
295 static const G4double prec = 1.e-03 * CLHEP::eV;
296 G4double coef, oldsum(0.), newsum(0.);
297 fMatNbOfIntervals = 0;
298
299 for (interval = 0; interval < interval2; ++interval) {
300 Emin = (*(*fMatSandiaMatrix)[fMatNbOfIntervals])[0] = tmp2[interval];
301
302 for (G4int k = 1; k < 5; ++k) {
303 (*(*fMatSandiaMatrix)[fMatNbOfIntervals])[k] = 0.;
304 }
305 newsum = 0.;
306
307 for (elm = 0; elm < NbElm; elm++) {
308 GetSandiaCofPerAtom(Z[elm], Emin + prec, fSandiaCofPerAtom);
309
310 for (G4int j = 1; j < 5; ++j) {
311 coef = NbOfAtomsPerVolume[elm] * fSandiaCofPerAtom[j - 1];
312 (*(*fMatSandiaMatrix)[fMatNbOfIntervals])[j] += coef;
313 newsum += std::abs(coef);
314 }
315 }
316 // check for null or redondant intervals
317
318 if (newsum != oldsum) {
319 oldsum = newsum;
320 ++fMatNbOfIntervals;
321 }
322 }
323 delete[] Z;
324 delete[] tmp1;
325 delete[] tmp2;
326
327 if (fVerbose > 0) {
328 G4cout << "G4SandiaTable::ComputeMatSandiaMatrix(), mat = " << fMaterial->GetName() << G4endl;
329
330 for (G4int i = 0; i < fMatNbOfIntervals; ++i) {
331 G4cout << i << "\t" << GetSandiaCofForMaterial(i, 0) / keV << " keV \t"
332 << GetSandiaCofForMaterial(i, 1) << "\t" << GetSandiaCofForMaterial(i, 2) << "\t"
333 << GetSandiaCofForMaterial(i, 3) << "\t" << GetSandiaCofForMaterial(i, 4) << G4endl;
334 }
335 }
336 }
337
338 ////////////////////////////////////////////////////////////////////////////////
339 //
340 // Sandia matrix for PAI models based on vectors ...
341
342 void G4SandiaTable::ComputeMatSandiaMatrixPAI()
343 {
344 G4int MaxIntervals = 0;
345 G4int elm, c, i, j, jj, k, k1, k2, c1, n1, z;
346
347 const auto noElm = (G4int)fMaterial->GetNumberOfElements();
348 const G4ElementVector* ElementVector = fMaterial->GetElementVector();
349
350 std::vector<G4int> Z(noElm); // Atomic number
351
352 for (elm = 0; elm < noElm; ++elm) {
353 z = G4lrint((*ElementVector)[elm]->GetZ());
354 if (z < 1) {
355 z = 1;
356 }
357 else if (z > 100) {
358 z = 100;
359 }
360 Z[elm] = z;
361 MaxIntervals += fNbOfIntervals[Z[elm]];
362 }
363 fMaxInterval = MaxIntervals + 2;
364
365 if (fVerbose > 0) {
366 G4cout << "G4SandiaTable::ComputeMatSandiaMatrixPAI: fMaxInterval = " << fMaxInterval << G4endl;
367 }
368
369 G4DataVector fPhotoAbsorptionCof0(fMaxInterval);
370 G4DataVector fPhotoAbsorptionCof1(fMaxInterval);
371 G4DataVector fPhotoAbsorptionCof2(fMaxInterval);
372 G4DataVector fPhotoAbsorptionCof3(fMaxInterval);
373 G4DataVector fPhotoAbsorptionCof4(fMaxInterval);
374
375 for (c = 0; c < fMaxInterval; ++c) // just in case
376 {
377 fPhotoAbsorptionCof0[c] = 0.;
378 fPhotoAbsorptionCof1[c] = 0.;
379 fPhotoAbsorptionCof2[c] = 0.;
380 fPhotoAbsorptionCof3[c] = 0.;
381 fPhotoAbsorptionCof4[c] = 0.;
382 }
383 c = 1;
384
385 for (i = 0; i < noElm; ++i) {
386 G4double I1 = fIonizationPotentials[Z[i]] * CLHEP::keV; // I1 in keV
387 n1 = 1;
388
389 for (j = 1; j < Z[i]; ++j) {
390 n1 += fNbOfIntervals[j];
391 }
392
393 G4int n2 = n1 + fNbOfIntervals[Z[i]];
394
395 for (k1 = n1; k1 < n2; ++k1) {
396 if (I1 > fSandiaTable[k1][0]) {
397 continue; // no ionization for energies smaller than I1 (first
398 } // ionisation potential)
399 break;
400 }
401 G4int flag = 0;
402
403 for (c1 = 1; c1 < c; ++c1) {
404 if (fPhotoAbsorptionCof0[c1] == I1) // this value already has existed
405 {
406 flag = 1;
407 break;
408 }
409 }
410 if (flag == 0) {
411 fPhotoAbsorptionCof0[c] = I1;
412 ++c;
413 }
414 for (k2 = k1; k2 < n2; ++k2) {
415 flag = 0;
416
417 for (c1 = 1; c1 < c; ++c1) {
418 if (fPhotoAbsorptionCof0[c1] == fSandiaTable[k2][0]) {
419 flag = 1;
420 break;
421 }
422 }
423 if (flag == 0) {
424 fPhotoAbsorptionCof0[c] = fSandiaTable[k2][0];
425 ++c;
426 }
427 }
428 } // end for(i)
429 // sort out
430
431 for (i = 1; i < c; ++i) {
432 for (j = i + 1; j < c; ++j) {
433 if (fPhotoAbsorptionCof0[i] > fPhotoAbsorptionCof0[j]) {
434 G4double tmp = fPhotoAbsorptionCof0[i];
435 fPhotoAbsorptionCof0[i] = fPhotoAbsorptionCof0[j];
436 fPhotoAbsorptionCof0[j] = tmp;
437 }
438 }
439 if (fVerbose > 0) {
440 G4cout << i << "\t energy = " << fPhotoAbsorptionCof0[i] << G4endl;
441 }
442 }
443 fMaxInterval = c;
444
445 const G4double* fractionW = fMaterial->GetFractionVector();
446
447 if (fVerbose > 0) {
448 for (i = 0; i < noElm; ++i) {
449 G4cout << i << " = elN, fraction = " << fractionW[i] << G4endl;
450 }
451 }
452
453 for (i = 0; i < noElm; ++i) {
454 n1 = 1;
455 G4double I1 = fIonizationPotentials[Z[i]] * keV;
456
457 for (j = 1; j < Z[i]; ++j) {
458 n1 += fNbOfIntervals[j];
459 }
460
461 G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1;
462
463 for (k = n1; k < n2; ++k) {
464 G4double B1 = fSandiaTable[k][0];
465 G4double B2 = fSandiaTable[k + 1][0];
466
467 for (G4int q = 1; q < fMaxInterval - 1; ++q) {
468 G4double E1 = fPhotoAbsorptionCof0[q];
469 G4double E2 = fPhotoAbsorptionCof0[q + 1];
470
471 if (fVerbose > 0) {
472 G4cout << "k = " << k << ", q = " << q << ", B1 = " << B1 << ", B2 = " << B2
473 << ", E1 = " << E1 << ", E2 = " << E2 << G4endl;
474 }
475 if (B1 > E1 || B2 < E2 || E1 < I1) {
476 if (fVerbose > 0) {
477 G4cout << "continue for: B1 = " << B1 << ", B2 = " << B2 << ", E1 = " << E1
478 << ", E2 = " << E2 << G4endl;
479 }
480 continue;
481 }
482 fPhotoAbsorptionCof1[q] += fSandiaTable[k][1] * fractionW[i];
483 fPhotoAbsorptionCof2[q] += fSandiaTable[k][2] * fractionW[i];
484 fPhotoAbsorptionCof3[q] += fSandiaTable[k][3] * fractionW[i];
485 fPhotoAbsorptionCof4[q] += fSandiaTable[k][4] * fractionW[i];
486 }
487 }
488 // Last interval
489
490 fPhotoAbsorptionCof1[fMaxInterval - 1] += fSandiaTable[k][1] * fractionW[i];
491 fPhotoAbsorptionCof2[fMaxInterval - 1] += fSandiaTable[k][2] * fractionW[i];
492 fPhotoAbsorptionCof3[fMaxInterval - 1] += fSandiaTable[k][3] * fractionW[i];
493 fPhotoAbsorptionCof4[fMaxInterval - 1] += fSandiaTable[k][4] * fractionW[i];
494 } // for(i)
495 c = 0; // Deleting of first intervals where all coefficients = 0
496
497 do {
498 ++c;
499
500 if (fPhotoAbsorptionCof1[c] != 0.0 || fPhotoAbsorptionCof2[c] != 0.0 ||
501 fPhotoAbsorptionCof3[c] != 0.0 || fPhotoAbsorptionCof4[c] != 0.0)
502 {
503 continue;
504 }
505
506 if (fVerbose > 0) {
507 G4cout << c << " = number with zero cofs" << G4endl;
508 }
509 for (jj = 2; jj < fMaxInterval; ++jj) {
510 fPhotoAbsorptionCof0[jj - 1] = fPhotoAbsorptionCof0[jj];
511 fPhotoAbsorptionCof1[jj - 1] = fPhotoAbsorptionCof1[jj];
512 fPhotoAbsorptionCof2[jj - 1] = fPhotoAbsorptionCof2[jj];
513 fPhotoAbsorptionCof3[jj - 1] = fPhotoAbsorptionCof3[jj];
514 fPhotoAbsorptionCof4[jj - 1] = fPhotoAbsorptionCof4[jj];
515 }
516 --fMaxInterval;
517 --c;
518 }
519 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
520 while (c < fMaxInterval - 1);
521
522 if (fPhotoAbsorptionCof0[fMaxInterval - 1] == 0.0) {
523 fMaxInterval--;
524 }
525
526 // create the sandia matrix for this material
527
528 fMatSandiaMatrixPAI = new G4OrderedTable();
529
530 G4double density = fMaterial->GetDensity();
531
532 for (i = 0; i < fMaxInterval; ++i) // -> G4units
533 {
534 fPhotoAbsorptionCof0[i + 1] *= funitc[0];
535 fPhotoAbsorptionCof1[i + 1] *= funitc[1] * density;
536 fPhotoAbsorptionCof2[i + 1] *= funitc[2] * density;
537 fPhotoAbsorptionCof3[i + 1] *= funitc[3] * density;
538 fPhotoAbsorptionCof4[i + 1] *= funitc[4] * density;
539 }
540 if (fLowerI1) {
541 if (fMaterial->GetName() == "G4_WATER") {
542 fMaxInterval += fH2OlowerInt;
543
544 for (i = 0; i < fMaxInterval; ++i) // init vector table
545 {
546 fMatSandiaMatrixPAI->push_back(new G4DataVector(5, 0.));
547 }
548 for (i = 0; i < fH2OlowerInt; ++i) {
549 (*(*fMatSandiaMatrixPAI)[i])[0] = fH2OlowerI1[i][0];
550 (*(*fMatSandiaMatrixPAI)[i])[1] = fH2OlowerI1[i][1];
551 (*(*fMatSandiaMatrixPAI)[i])[2] = fH2OlowerI1[i][2];
552 (*(*fMatSandiaMatrixPAI)[i])[3] = fH2OlowerI1[i][3];
553 (*(*fMatSandiaMatrixPAI)[i])[4] = fH2OlowerI1[i][4];
554 }
555 for (i = fH2OlowerInt; i < fMaxInterval; ++i) {
556 (*(*fMatSandiaMatrixPAI)[i])[0] = fPhotoAbsorptionCof0[i + 1 - fH2OlowerInt];
557 (*(*fMatSandiaMatrixPAI)[i])[1] = fPhotoAbsorptionCof1[i + 1 - fH2OlowerInt];
558 (*(*fMatSandiaMatrixPAI)[i])[2] = fPhotoAbsorptionCof2[i + 1 - fH2OlowerInt];
559 (*(*fMatSandiaMatrixPAI)[i])[3] = fPhotoAbsorptionCof3[i + 1 - fH2OlowerInt];
560 (*(*fMatSandiaMatrixPAI)[i])[4] = fPhotoAbsorptionCof4[i + 1 - fH2OlowerInt];
561 }
562 }
563 }
564 else {
565 for (i = 0; i < fMaxInterval; ++i) // init vector table
566 {
567 fMatSandiaMatrixPAI->push_back(new G4DataVector(5, 0.));
568 }
569 for (i = 0; i < fMaxInterval; ++i) {
570 (*(*fMatSandiaMatrixPAI)[i])[0] = fPhotoAbsorptionCof0[i + 1];
571 (*(*fMatSandiaMatrixPAI)[i])[1] = fPhotoAbsorptionCof1[i + 1]; // *density;
572 (*(*fMatSandiaMatrixPAI)[i])[2] = fPhotoAbsorptionCof2[i + 1]; // *density;
573 (*(*fMatSandiaMatrixPAI)[i])[3] = fPhotoAbsorptionCof3[i + 1]; // *density;
574 (*(*fMatSandiaMatrixPAI)[i])[4] = fPhotoAbsorptionCof4[i + 1]; // *density;
575 }
576 }
577 // --fMaxInterval;
578 // to avoid duplicate at 500 keV or extra zeros in last interval
579
580 if (fVerbose > 0) {
581 G4cout << "G4SandiaTable::ComputeMatSandiaMatrixPAI(), mat = " << fMaterial->GetName()
582 << G4endl;
583
584 for (i = 0; i < fMaxInterval; ++i) {
585 G4cout << i << "\t" << GetSandiaMatTablePAI(i, 0) / keV << " keV \t"
586 << GetSandiaMatTablePAI(i, 1) << "\t" << GetSandiaMatTablePAI(i, 2) << "\t"
587 << GetSandiaMatTablePAI(i, 3) << "\t" << GetSandiaMatTablePAI(i, 4) << G4endl;
588 }
589 }
590 return;
591 }
592
593 ////////////////////////////////////////////////////////////////////////////////
594 // Methods for PAI model only
595 //
596
597 G4SandiaTable::G4SandiaTable(G4int matIndex)
598 {
599 fMaterial = nullptr;
600 fMatNbOfIntervals = 0;
601 fMatSandiaMatrix = nullptr;
602 fMatSandiaMatrixPAI = nullptr;
603 fPhotoAbsorptionCof = nullptr;
604
605 fMaxInterval = 0;
606 fVerbose = 0;
607 fLowerI1 = false;
608
609 fSandiaCofPerAtom.resize(4, 0.0);
610
611 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
612 auto numberOfMat = (G4int)G4Material::GetNumberOfMaterials();
613
614 if (matIndex >= 0 && matIndex < numberOfMat) {
615 fMaterial = (*theMaterialTable)[matIndex];
616 }
617 else {
618 G4Exception(
619 "G4SandiaTable::G4SandiaTable(G4int matIndex)", "mat401", FatalException, "wrong matIndex");
620 }
621 }
622
623 ////////////////////////////////////////////////////////////////////////////////
624
625 G4SandiaTable::G4SandiaTable()
626 {
627 fMaterial = nullptr;
628 fMatNbOfIntervals = 0;
629 fMatSandiaMatrix = nullptr;
630 fMatSandiaMatrixPAI = nullptr;
631 fPhotoAbsorptionCof = nullptr;
632
633 fMaxInterval = 0;
634 fVerbose = 0;
635 fLowerI1 = false;
636
637 fSandiaCofPerAtom.resize(4, 0.0);
638 }
639
640 ////////////////////////////////////////////////////////////////////////////////
641
642 void G4SandiaTable::Initialize(const G4Material* mat)
643 {
644 fMaterial = mat;
645 ComputeMatSandiaMatrixPAI();
646 }
647
648 ////////////////////////////////////////////////////////////////////////////////
649
650 G4int G4SandiaTable::GetMaxInterval() const { return fMaxInterval; }
651
652 ////////////////////////////////////////////////////////////////////////////////
653
654 G4double** G4SandiaTable::GetPointerToCof()
655 {
656 if (fPhotoAbsorptionCof == nullptr) {
657 ComputeMatTable();
658 }
659 return fPhotoAbsorptionCof;
660 }
661
662 ////////////////////////////////////////////////////////////////////////////////
663
664 void G4SandiaTable::SandiaSwap(G4double** da, G4int i, G4int j)
665 {
666 G4double tmp = da[i][0];
667 da[i][0] = da[j][0];
668 da[j][0] = tmp;
669 }
670
671 ////////////////////////////////////////////////////////////////////////////////
672
673 G4double G4SandiaTable::GetPhotoAbsorpCof(G4int i, G4int j) const
674 {
675 return fPhotoAbsorptionCof[i][j] * funitc[j];
676 }
677
678 ////////////////////////////////////////////////////////////////////////////////
679 //
680 // Bubble sorting of left energy interval in SandiaTable in ascening order
681 //
682
683 void G4SandiaTable::SandiaSort(G4double** da, G4int sz)
684 {
685 for (G4int i = 1; i < sz; ++i) {
686 for (G4int j = i + 1; j < sz; ++j) {
687 if (da[i][0] > da[j][0]) {
688 SandiaSwap(da, i, j);
689 }
690 }
691 }
692 }
693
694 ////////////////////////////////////////////////////////////////////////////////
695 //
696 // SandiaIntervals
697 //
698
699 G4int G4SandiaTable::SandiaIntervals(G4int Z[], G4int el)
700 {
701 G4int c, i, flag = 0, n1 = 1;
702 G4int j, c1, k1, k2;
703 G4double I1;
704 fMaxInterval = 0;
705
706 for (i = 0; i < el; ++i) {
707 fMaxInterval += fNbOfIntervals[Z[i]];
708 }
709
710 fMaxInterval += 2;
711
712 if (fVerbose > 0) {
713 G4cout << "begin sanInt, fMaxInterval = " << fMaxInterval << G4endl;
714 }
715
716 fPhotoAbsorptionCof = new G4double*[fMaxInterval];
717
718 for (i = 0; i < fMaxInterval; ++i) {
719 fPhotoAbsorptionCof[i] = new G4double[5];
720 }
721 // for(c = 0; c < fIntervalLimit; ++c) // just in case
722
723 for (c = 0; c < fMaxInterval; ++c) {
724 fPhotoAbsorptionCof[c][0] = 0.;
725 }
726
727 c = 1;
728
729 for (i = 0; i < el; ++i) {
730 I1 = fIonizationPotentials[Z[i]] * keV; // First ionization
731 n1 = 1; // potential in keV
732
733 for (j = 1; j < Z[i]; ++j) {
734 n1 += fNbOfIntervals[j];
735 }
736
737 G4int n2 = n1 + fNbOfIntervals[Z[i]];
738
739 for (k1 = n1; k1 < n2; k1++) {
740 if (I1 > fSandiaTable[k1][0]) {
741 continue; // no ionization for energies smaller than I1 (first
742 } // ionisation potential)
743 break;
744 }
745 flag = 0;
746
747 for (c1 = 1; c1 < c; c1++) {
748 if (fPhotoAbsorptionCof[c1][0] == I1) // this value already has existed
749 {
750 flag = 1;
751 break;
752 }
753 }
754 if (flag == 0) {
755 fPhotoAbsorptionCof[c][0] = I1;
756 ++c;
757 }
758 for (k2 = k1; k2 < n2; k2++) {
759 flag = 0;
760
761 for (c1 = 1; c1 < c; c1++) {
762 if (fPhotoAbsorptionCof[c1][0] == fSandiaTable[k2][0]) {
763 flag = 1;
764 break;
765 }
766 }
767 if (flag == 0) {
768 fPhotoAbsorptionCof[c][0] = fSandiaTable[k2][0];
769 if (fVerbose > 0) {
770 G4cout << "sanInt, c = " << c << ", E_c = " << fPhotoAbsorptionCof[c][0] << G4endl;
771 }
772 ++c;
773 }
774 }
775 } // end for(i)
776
777 SandiaSort(fPhotoAbsorptionCof, c);
778 fMaxInterval = c;
779 if (fVerbose > 0) {
780 G4cout << "end SanInt, fMaxInterval = " << fMaxInterval << G4endl;
781 }
782 return c;
783 }
784
785 ////////////////////////////////////////////////////////////////////////////..
786 //
787 // SandiaMixing
788 //
789
790 G4int G4SandiaTable::SandiaMixing(G4int Z[], const G4double fractionW[], G4int el, G4int mi)
791 {
792 G4int i, j, n1, k, c = 1, jj, kk;
793 G4double I1, B1, B2, E1, E2;
794
795 for (i = 0; i < mi; ++i) {
796 for (j = 1; j < 5; ++j) {
797 fPhotoAbsorptionCof[i][j] = 0.;
798 }
799 }
800 for (i = 0; i < el; ++i) {
801 n1 = 1;
802 I1 = fIonizationPotentials[Z[i]] * keV;
803
804 for (j = 1; j < Z[i]; ++j) {
805 n1 += fNbOfIntervals[j];
806 }
807
808 G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1;
809
810 for (k = n1; k < n2; ++k) {
811 B1 = fSandiaTable[k][0];
812 B2 = fSandiaTable[k + 1][0];
813
814 for (c = 1; c < mi - 1; ++c) {
815 E1 = fPhotoAbsorptionCof[c][0];
816 E2 = fPhotoAbsorptionCof[c + 1][0];
817
818 if (B1 > E1 || B2 < E2 || E1 < I1) {
819 continue;
820 }
821
822 for (j = 1; j < 5; ++j) {
823 fPhotoAbsorptionCof[c][j] += fSandiaTable[k][j] * fractionW[i];
824 if (fVerbose > 0) {
825 G4cout << "c=" << c << "; j=" << j << "; fST=" << fSandiaTable[k][j]
826 << "; frW=" << fractionW[i] << G4endl;
827 }
828 }
829 }
830 }
831 for (j = 1; j < 5; ++j) // Last interval
832 {
833 fPhotoAbsorptionCof[mi - 1][j] += fSandiaTable[k][j] * fractionW[i];
834 if (fVerbose > 0) {
835 G4cout << "mi-1=" << mi - 1 << "; j=" << j << "; fST=" << fSandiaTable[k][j]
836 << "; frW=" << fractionW[i] << G4endl;
837 }
838 }
839 } // for(i)
840 c = 0; // Deleting of first intervals where all coefficients = 0
841
842 do {
843 ++c;
844
845 if (fPhotoAbsorptionCof[c][1] != 0.0 || fPhotoAbsorptionCof[c][2] != 0.0 ||
846 fPhotoAbsorptionCof[c][3] != 0.0 || fPhotoAbsorptionCof[c][4] != 0.0)
847 {
848 continue;
849 }
850
851 for (jj = 2; jj < mi; ++jj) {
852 for (kk = 0; kk < 5; ++kk) {
853 fPhotoAbsorptionCof[jj - 1][kk] = fPhotoAbsorptionCof[jj][kk];
854 }
855 }
856 mi--;
857 c--;
858 }
859 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
860 while (c < mi - 1);
861
862 if (fVerbose > 0) {
863 G4cout << "end SanMix, mi = " << mi << G4endl;
864 }
865
866 return mi;
867 }
868
869 ////////////////////////////////////////////////////////////////////////////////
870
871 G4int G4SandiaTable::GetMatNbOfIntervals() const { return fMatNbOfIntervals; }
872
873 ////////////////////////////////////////////////////////////////////////////////
874
875 G4double G4SandiaTable::GetSandiaPerAtom(G4int Z, G4int interval, G4int j) const
876 {
877 #ifdef G4VERBOSE
878 if (Z < 1 || Z > 100) {
879 Z = PrintErrorZ(Z, "GetSandiaPerAtom");
880 }
881 if (interval < 0 || interval >= fNbOfIntervals[Z]) {
882 PrintErrorV("GetSandiaPerAtom");
883 interval = (interval < 0) ? 0 : fNbOfIntervals[Z] - 1;
884 }
885 if (j < 0 || j > 4) {
886 PrintErrorV("GetSandiaPerAtom");
887 j = (j < 0) ? 0 : 4;
888 }
889 #endif
890 G4int row = fCumulInterval[Z - 1] + interval;
891 G4double x = fSandiaTable[row][0] * CLHEP::keV;
892 if (j > 0) {
893 x = Z * CLHEP::amu / fZtoAratio[Z] * fSandiaTable[row][j] * funitc[j];
894 }
895 return x;
896 }
897
898 ////////////////////////////////////////////////////////////////////////////////
899
900 G4double G4SandiaTable::GetSandiaCofForMaterial(G4int interval, G4int j) const
901 {
902 #ifdef G4VERBOSE
903 if (interval < 0 || interval >= fMatNbOfIntervals) {
904 PrintErrorV("GetSandiaCofForMaterial");
905 interval = (interval < 0) ? 0 : fMatNbOfIntervals - 1;
906 }
907 if (j < 0 || j > 4) {
908 PrintErrorV("GetSandiaCofForMaterial");
909 j = (j < 0) ? 0 : 4;
910 }
911 #endif
912 return ((*(*fMatSandiaMatrix)[interval])[j]);
913 }
914
915 ////////////////////////////////////////////////////////////////////////////////
916
917 const G4double* G4SandiaTable::GetSandiaCofForMaterial(G4double energy) const
918 {
919 G4int interval = 0;
920 if (energy > (*(*fMatSandiaMatrix)[0])[0]) {
921 interval = fMatNbOfIntervals - 1;
922 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
923 while ((interval > 0) && (energy < (*(*fMatSandiaMatrix)[interval])[0])) {
924 --interval;
925 }
926 }
927 return &((*(*fMatSandiaMatrix)[interval])[1]);
928 }
929
930 ////////////////////////////////////////////////////////////////////////////////
931
932 G4double G4SandiaTable::GetSandiaMatTable(G4int interval, G4int j) const
933 {
934 #ifdef G4VERBOSE
935 if (interval < 0 || interval >= fMatNbOfIntervals) {
936 PrintErrorV("GetSandiaCofForMaterial");
937 interval = (interval < 0) ? 0 : fMatNbOfIntervals - 1;
938 }
939 if (j < 0 || j > 4) {
940 PrintErrorV("GetSandiaCofForMaterial");
941 j = (j < 0) ? 0 : 4;
942 }
943 #endif
944 return ((*(*fMatSandiaMatrix)[interval])[j]) * funitc[j];
945 }
946
947 ////////////////////////////////////////////////////////////////////////////////
948
949 G4double G4SandiaTable::GetSandiaMatTablePAI(G4int interval, G4int j) const
950 {
951 #ifdef G4VERBOSE
952 if (interval < 0 || interval >= fMaxInterval) {
953 PrintErrorV("GetSandiaCofForMaterialPAI");
954 interval = (interval < 0) ? 0 : fMaxInterval - 1;
955 }
956 if (j < 0 || j > 4) {
957 PrintErrorV("GetSandiaCofForMaterialPAI");
958 j = (j < 0) ? 0 : 4;
959 }
960 #endif
961 return ((*(*fMatSandiaMatrixPAI)[interval])[j]);
962 }
963
964 ////////////////////////////////////////////////////////////////////////////////
965 //
966 // Sandia interval and mixing calculations for materialCutsCouple constructor
967 //
968
969 void G4SandiaTable::ComputeMatTable()
970 {
971 G4int elm, c, i, j, jj, k, kk, k1, k2, c1, n1;
972
973 const auto noElm = (G4int)fMaterial->GetNumberOfElements();
974 const G4ElementVector* ElementVector = fMaterial->GetElementVector();
975 auto Z = new G4int[noElm]; // Atomic number
976
977 fMaxInterval = 0;
978 for (elm = 0; elm < noElm; ++elm) {
979 Z[elm] = (*ElementVector)[elm]->GetZasInt();
980 fMaxInterval += fNbOfIntervals[Z[elm]];
981 }
982 fMaxInterval += 2;
983
984 // G4cout<<"fMaxInterval = "<<fMaxInterval<<G4endl;
985
986 fPhotoAbsorptionCof = new G4double*[fMaxInterval];
987
988 for (i = 0; i < fMaxInterval; ++i) {
989 fPhotoAbsorptionCof[i] = new G4double[5];
990 }
991
992 // for(c = 0; c < fIntervalLimit; ++c) // just in case
993
994 for (c = 0; c < fMaxInterval; ++c) // just in case
995 {
996 fPhotoAbsorptionCof[c][0] = 0.;
997 }
998 c = 1;
999
1000 for (i = 0; i < noElm; ++i) {
1001 G4double I1 = fIonizationPotentials[Z[i]] * keV; // First ionization
1002 n1 = 1; // potential in keV
1003
1004 for (j = 1; j < Z[i]; ++j) {
1005 n1 += fNbOfIntervals[j];
1006 }
1007 G4int n2 = n1 + fNbOfIntervals[Z[i]];
1008
1009 for (k1 = n1; k1 < n2; ++k1) {
1010 if (I1 > fSandiaTable[k1][0]) {
1011 continue; // no ionization for energies smaller than I1 (first
1012 } // ionisation potential)
1013 break;
1014 }
1015 G4int flag = 0;
1016
1017 for (c1 = 1; c1 < c; ++c1) {
1018 if (fPhotoAbsorptionCof[c1][0] == I1) // this value already has existed
1019 {
1020 flag = 1;
1021 break;
1022 }
1023 }
1024 if (flag == 0) {
1025 fPhotoAbsorptionCof[c][0] = I1;
1026 ++c;
1027 }
1028 for (k2 = k1; k2 < n2; ++k2) {
1029 flag = 0;
1030
1031 for (c1 = 1; c1 < c; ++c1) {
1032 if (fPhotoAbsorptionCof[c1][0] == fSandiaTable[k2][0]) {
1033 flag = 1;
1034 break;
1035 }
1036 }
1037 if (flag == 0) {
1038 fPhotoAbsorptionCof[c][0] = fSandiaTable[k2][0];
1039 ++c;
1040 }
1041 }
1042 } // end for(i)
1043
1044 SandiaSort(fPhotoAbsorptionCof, c);
1045 fMaxInterval = c;
1046
1047 const G4double* fractionW = fMaterial->GetFractionVector();
1048
1049 for (i = 0; i < fMaxInterval; ++i) {
1050 for (j = 1; j < 5; ++j) {
1051 fPhotoAbsorptionCof[i][j] = 0.;
1052 }
1053 }
1054 for (i = 0; i < noElm; ++i) {
1055 n1 = 1;
1056 G4double I1 = fIonizationPotentials[Z[i]] * keV;
1057
1058 for (j = 1; j < Z[i]; ++j) {
1059 n1 += fNbOfIntervals[j];
1060 }
1061 G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1;
1062
1063 for (k = n1; k < n2; ++k) {
1064 G4double B1 = fSandiaTable[k][0];
1065 G4double B2 = fSandiaTable[k + 1][0];
1066 for (G4int q = 1; q < fMaxInterval - 1; q++) {
1067 G4double E1 = fPhotoAbsorptionCof[q][0];
1068 G4double E2 = fPhotoAbsorptionCof[q + 1][0];
1069 if (B1 > E1 || B2 < E2 || E1 < I1) {
1070 continue;
1071 }
1072 for (j = 1; j < 5; ++j) {
1073 fPhotoAbsorptionCof[q][j] += fSandiaTable[k][j] * fractionW[i];
1074 }
1075 }
1076 }
1077 for (j = 1; j < 5; ++j) // Last interval
1078 {
1079 fPhotoAbsorptionCof[fMaxInterval - 1][j] += fSandiaTable[k][j] * fractionW[i];
1080 }
1081 } // for(i)
1082
1083 c = 0; // Deleting of first intervals where all coefficients = 0
1084
1085 do {
1086 ++c;
1087
1088 if (fPhotoAbsorptionCof[c][1] != 0.0 || fPhotoAbsorptionCof[c][2] != 0.0 ||
1089 fPhotoAbsorptionCof[c][3] != 0.0 || fPhotoAbsorptionCof[c][4] != 0.0)
1090 {
1091 continue;
1092 }
1093
1094 for (jj = 2; jj < fMaxInterval; ++jj) {
1095 for (kk = 0; kk < 5; ++kk) {
1096 fPhotoAbsorptionCof[jj - 1][kk] = fPhotoAbsorptionCof[jj][kk];
1097 }
1098 }
1099 --fMaxInterval;
1100 --c;
1101 }
1102 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
1103 while (c < fMaxInterval - 1);
1104
1105 // create the sandia matrix for this material
1106
1107 --fMaxInterval; // vmg 20.11.10
1108
1109 fMatSandiaMatrix = new G4OrderedTable();
1110
1111 for (i = 0; i < fMaxInterval; ++i) {
1112 fMatSandiaMatrix->push_back(new G4DataVector(5, 0.));
1113 }
1114 for (i = 0; i < fMaxInterval; ++i) {
1115 for (j = 0; j < 5; ++j) {
1116 (*(*fMatSandiaMatrix)[i])[j] = fPhotoAbsorptionCof[i + 1][j];
1117 }
1118 }
1119 fMatNbOfIntervals = fMaxInterval;
1120
1121 if (fVerbose > 0) {
1122 G4cout << "vmg, G4SandiaTable::ComputeMatTable(), mat = " << fMaterial->GetName() << G4endl;
1123
1124 for (i = 0; i < fMaxInterval; ++i) {
1125 G4cout << i << "\t" << GetSandiaCofForMaterial(i, 0) / keV << " keV \t"
1126 << this->GetSandiaCofForMaterial(i, 1) << "\t" << this->GetSandiaCofForMaterial(i, 2)
1127 << "\t" << this->GetSandiaCofForMaterial(i, 3) << "\t"
1128 << this->GetSandiaCofForMaterial(i, 4) << G4endl;
1129 }
1130 }
1131 delete[] Z;
1132 return;
1133 }
1134