102 : fMaterial(nullptr),fMatSandiaMatrix(nullptr),
103 fMatSandiaMatrixPAI(nullptr),fPhotoAbsorptionCof(nullptr)
136 std::vector<G4double>& coeff)
const
139 if(Z < 1 || Z > 100) {
140 Z = PrintErrorZ(Z,
"GetSandiaCofPerAtom");
142 if(4 > coeff.size()) {
143 PrintErrorV(
"GetSandiaCofPerAtom(): input vector is resized");
152 if (energy <= Emin) {
159 while ((interval>0) && (energy<
fSandiaTable[row][0]*CLHEP::keV)) {
168 coeff[1]=AoverAvo*
funitc[2]*fSandiaTable[row][2];
169 coeff[2]=AoverAvo*
funitc[3]*fSandiaTable[row][3];
170 coeff[3]=AoverAvo*
funitc[4]*fSandiaTable[row][4];
177 std::vector<G4double>& coeff)
const
180 if(4 > coeff.size()) {
181 PrintErrorV(
"GetSandiaCofWater: input vector is resized");
189 if(energy >=
fH2OlowerI1[i][0]*CLHEP::keV) {
break; }
193 coeff[1]=
funitc[2]*fH2OlowerI1[i][2];
194 coeff[2]=
funitc[3]*fH2OlowerI1[i][3];
195 coeff[3]=
funitc[4]*fH2OlowerI1[i][4];
217 if(Z < 1 || Z > 100) {
218 Z = PrintErrorZ(Z,
"GetSandiaCofPerAtom");
232 ed <<
"Atomic number out of range Z= " << Z <<
"; closest value is used";
234 return (Z > 100) ? 100 : 1;
239 void G4SandiaTable::PrintErrorV(
const G4String& ss)
258 G4int MaxIntervals = 0;
263 for (elm = 0; elm < NbElm; ++elm)
265 z =
G4lrint((*ElementVector)[elm]->GetZ());
267 else if(z > 100) { z = 100; }
278 for (elm = 0; elm < NbElm; ++elm)
300 for (
G4int i1 = 0; i1 < MaxIntervals; ++i1)
305 tmp2[interval2] =
Emin;
309 for (
G4int j1 = 0; j1 < MaxIntervals; ++j1)
311 if (tmp1[j1] <= Emin) { tmp1[j1] =
DBL_MAX; }
321 for (interval = 0; interval < interval2; ++interval)
331 G4double coef, oldsum(0.), newsum(0.);
334 for ( interval = 0; interval < interval2; ++interval)
338 for (
G4int k = 1; k < 5; ++k ) {
343 for ( elm = 0; elm < NbElm; elm++ )
347 for (
G4int j = 1; j < 5; ++j )
351 newsum += std::abs(coef);
364 G4cout<<
"G4SandiaTable::ComputeMatSandiaMatrix(), mat = "
384 G4int MaxIntervals = 0;
385 G4int elm, c, i, j, jj, k, k1, k2,
c1, n1,
z;
390 std::vector<G4int>
Z(noElm);
392 for ( elm = 0; elm < noElm; elm++ )
394 z =
G4lrint((*ElementVector)[elm]->GetZ());
396 else if(z > 100) { z = 100; }
404 G4cout<<
"G4SandiaTable::ComputeMatSandiaMatrixPAI: fMaxInterval = "
416 fPhotoAbsorptionCof0[c] = 0.;
417 fPhotoAbsorptionCof1[c] = 0.;
418 fPhotoAbsorptionCof2[c] = 0.;
419 fPhotoAbsorptionCof3[c] = 0.;
420 fPhotoAbsorptionCof4[c] = 0.;
424 for(i = 0; i < noElm; ++i)
433 for( k1 = n1; k1 < n2; k1++ )
443 for( c1 = 1; c1 < c; c1++ )
445 if( fPhotoAbsorptionCof0[c1] == I1 )
453 fPhotoAbsorptionCof0[c] = I1;
456 for( k2 = k1; k2 < n2; k2++ )
460 for( c1 = 1; c1 < c; c1++ )
477 for( i = 1; i < c; ++i )
479 for( j = i + 1; j < c; ++j )
481 if( fPhotoAbsorptionCof0[i] > fPhotoAbsorptionCof0[j] )
484 fPhotoAbsorptionCof0[i] = fPhotoAbsorptionCof0[j];
485 fPhotoAbsorptionCof0[j] =
tmp;
490 G4cout<<i<<
"\t energy = "<<fPhotoAbsorptionCof0[i]<<
G4endl;
499 for( i = 0; i < noElm; ++i )
503 for( i = 0; i < noElm; ++i )
512 for(k = n1; k < n2; ++k)
517 for(
G4int q = 1; q < fMaxInterval-1; q++)
519 G4double E1 = fPhotoAbsorptionCof0[q];
520 G4double E2 = fPhotoAbsorptionCof0[q+1];
524 G4cout<<
"k = "<<k<<
", q = "<<q<<
", B1 = "<<B1<<
", B2 = "<<B2
525 <<
", E1 = "<<E1<<
", E2 = "<<E2<<
G4endl;
527 if( B1 > E1 || B2 < E2 || E1 < I1 )
531 G4cout<<
"continue for: B1 = "<<B1<<
", B2 = "<<B2<<
", E1 = "
532 <<E1<<
", E2 = "<<E2<<
G4endl;
536 fPhotoAbsorptionCof1[q] +=
fSandiaTable[k][1]*fractionW[i];
537 fPhotoAbsorptionCof2[q] +=
fSandiaTable[k][2]*fractionW[i];
538 fPhotoAbsorptionCof3[q] +=
fSandiaTable[k][3]*fractionW[i];
539 fPhotoAbsorptionCof4[q] +=
fSandiaTable[k][4]*fractionW[i];
544 fPhotoAbsorptionCof1[fMaxInterval-1] +=
fSandiaTable[k][1]*fractionW[i];
545 fPhotoAbsorptionCof2[fMaxInterval-1] +=
fSandiaTable[k][2]*fractionW[i];
546 fPhotoAbsorptionCof3[fMaxInterval-1] +=
fSandiaTable[k][3]*fractionW[i];
547 fPhotoAbsorptionCof4[fMaxInterval-1] +=
fSandiaTable[k][4]*fractionW[i];
555 if( fPhotoAbsorptionCof1[c] != 0.0 ||
556 fPhotoAbsorptionCof2[c] != 0.0 ||
557 fPhotoAbsorptionCof3[c] != 0.0 ||
558 fPhotoAbsorptionCof4[c] != 0.0 )
continue;
566 fPhotoAbsorptionCof0[jj-1] = fPhotoAbsorptionCof0[jj];
567 fPhotoAbsorptionCof1[jj-1] = fPhotoAbsorptionCof1[jj];
568 fPhotoAbsorptionCof2[jj-1] = fPhotoAbsorptionCof2[jj];
569 fPhotoAbsorptionCof3[jj-1] = fPhotoAbsorptionCof3[jj];
570 fPhotoAbsorptionCof4[jj-1] = fPhotoAbsorptionCof4[jj];
576 while( c < fMaxInterval - 1 );
578 if( fPhotoAbsorptionCof0[fMaxInterval-1] == 0.0 ) fMaxInterval--;
588 fPhotoAbsorptionCof0[i+1] *=
funitc[0];
589 fPhotoAbsorptionCof1[i+1] *=
funitc[1]*density;
590 fPhotoAbsorptionCof2[i+1] *=
funitc[2]*density;
591 fPhotoAbsorptionCof3[i+1] *=
funitc[3]*density;
592 fPhotoAbsorptionCof4[i+1] *=
funitc[4]*density;
606 (*(*fMatSandiaMatrixPAI)[i])[0] =
fH2OlowerI1[i][0];
607 (*(*fMatSandiaMatrixPAI)[i])[1] =
fH2OlowerI1[i][1];
608 (*(*fMatSandiaMatrixPAI)[i])[2] =
fH2OlowerI1[i][2];
609 (*(*fMatSandiaMatrixPAI)[i])[3] =
fH2OlowerI1[i][3];
610 (*(*fMatSandiaMatrixPAI)[i])[4] =
fH2OlowerI1[i][4];
614 (*(*fMatSandiaMatrixPAI)[i])[0] = fPhotoAbsorptionCof0[i+1-fH2OlowerInt];
615 (*(*fMatSandiaMatrixPAI)[i])[1] = fPhotoAbsorptionCof1[i+1-fH2OlowerInt];
616 (*(*fMatSandiaMatrixPAI)[i])[2] = fPhotoAbsorptionCof2[i+1-fH2OlowerInt];
617 (*(*fMatSandiaMatrixPAI)[i])[3] = fPhotoAbsorptionCof3[i+1-fH2OlowerInt];
618 (*(*fMatSandiaMatrixPAI)[i])[4] = fPhotoAbsorptionCof4[i+1-fH2OlowerInt];
630 (*(*fMatSandiaMatrixPAI)[i])[0] = fPhotoAbsorptionCof0[i+1];
631 (*(*fMatSandiaMatrixPAI)[i])[1] = fPhotoAbsorptionCof1[i+1];
632 (*(*fMatSandiaMatrixPAI)[i])[2] = fPhotoAbsorptionCof2[i+1];
633 (*(*fMatSandiaMatrixPAI)[i])[3] = fPhotoAbsorptionCof3[i+1];
634 (*(*fMatSandiaMatrixPAI)[i])[4] = fPhotoAbsorptionCof4[i+1];
642 G4cout<<
"G4SandiaTable::ComputeMatSandiaMatrixPAI(), mat = "
678 if ( matIndex >= 0 && matIndex < numberOfMat)
680 fMaterial = (*theMaterialTable)[matIndex];
684 G4Exception(
"G4SandiaTable::G4SandiaTable(G4int matIndex)",
"mat401",
736 da[i][0] = da[j][0] ;
755 for(
G4int i = 1;i < sz; ++i )
757 for(
G4int j = i + 1;j < sz; ++j )
771 G4int c, i, flag = 0, n1 = 1;
795 for( i = 0; i < el; ++i )
804 for( k1 = n1; k1 < n2; k1++ )
814 for( c1 = 1; c1 < c; c1++ )
827 for( k2 = k1; k2 < n2; k2++ )
831 for( c1 = 1; c1 < c; c1++ )
854 G4cout<<
"end SanInt, fMaxInterval = "<<fMaxInterval<<
G4endl;
870 G4int i, j, n1, k, c=1, jj, kk;
873 for( i = 0; i < mi; ++i )
877 for( i = 0; i < el; ++i )
886 for( k = n1; k < n2; ++k )
891 for( c = 1; c < mi-1; ++c )
896 if( B1 > E1 || B2 < E2 || E1 < I1 )
continue;
898 for( j = 1; j < 5; ++j )
904 <<
"; frW="<<fractionW[i]<<
G4endl;
909 for( j = 1; j < 5; ++j )
915 <<
"; frW="<<fractionW[i]<<
G4endl;
930 for( jj = 2; jj < mi; ++jj )
932 for( kk = 0; kk < 5; ++kk ) {
960 if(Z < 1 || Z > 100) {
961 Z = PrintErrorZ(Z,
"GetSandiaPerAtom");
964 PrintErrorV(
"GetSandiaPerAtom");
968 PrintErrorV(
"GetSandiaPerAtom");
987 PrintErrorV(
"GetSandiaCofForMaterial");
991 PrintErrorV(
"GetSandiaCofForMaterial");
1010 return &((*(*fMatSandiaMatrix)[interval])[1]);
1020 PrintErrorV(
"GetSandiaCofForMaterial");
1024 PrintErrorV(
"GetSandiaCofForMaterial");
1038 PrintErrorV(
"GetSandiaCofForMaterialPAI");
1042 PrintErrorV(
"GetSandiaCofForMaterialPAI");
1056 G4int MaxIntervals = 0;
1057 G4int elm, c, i, j, jj, k, kk, k1, k2,
c1, n1;
1063 for (elm = 0; elm<noElm; ++elm)
1065 Z[elm] = (
G4int)(*ElementVector)[elm]->GetZ();
1091 for(i = 0; i < noElm; ++i)
1096 for(j = 1; j < Z[i]; ++j)
1102 for(k1 = n1; k1 < n2; ++k1)
1112 for(c1 = 1; c1 < c; ++
c1)
1125 for(k2 = k1; k2 < n2; ++k2)
1129 for(c1 = 1; c1 < c; ++
c1)
1154 for(i = 0; i < noElm; ++i)
1159 for(j = 1; j < Z[i]; ++j)
1165 for(k = n1; k < n2; ++k)
1169 for(
G4int q = 1; q < fMaxInterval-1; q++)
1173 if(B1 > E1 || B2 < E2 || E1 < I1)
1177 for(j = 1; j < 5; ++j)
1183 for(j = 1; j < 5; ++j)
1203 for(kk = 0; kk < 5; ++kk)
1212 while( c < fMaxInterval - 1 );
1226 for( j = 0; j < 5; ++j )
1235 G4cout<<
"vmg, G4SandiaTable::ComputeMatTable(), mat = "
G4int SandiaMixing(G4int Z[], const G4double *fractionW, G4int el, G4int mi)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static const G4double Emin
const G4double * GetFractionVector() const
std::ostringstream G4ExceptionDescription
G4double ** GetPointerToCof()
static size_t GetNumberOfMaterials()
static G4MaterialTable * GetMaterialTable()
static constexpr double keV
static constexpr double g
const G4ElementVector * GetElementVector() const
static const G4double fIonizationPotentials[101]
std::vector< G4double > fSandiaCofPerAtom
void SandiaSwap(G4double **da, G4int i, G4int j)
void ComputeMatSandiaMatrix()
const G4double * GetVecNbOfAtomsPerVolume() const
G4double GetSandiaPerAtom(G4int Z, G4int, G4int) const
static const G4int fNbOfIntervals[101]
G4int GetMaxInterval() const
void SandiaSort(G4double **da, G4int sz)
G4int GetMatNbOfIntervals() const
const G4String & GetName() const
static const char sss[MAX_N_PAR+2]
void GetSandiaCofWater(G4double energy, std::vector< G4double > &coeff) const
static G4double GetZtoA(G4int Z)
void Initialize(G4Material *)
G4double ** fPhotoAbsorptionCof
void GetSandiaCofPerAtom(G4int Z, G4double energy, std::vector< G4double > &coeff) const
G4int SandiaIntervals(G4int Z[], G4int el)
static constexpr double eV
static const G4int fH2OlowerInt
G4double GetSandiaMatTablePAI(G4int, G4int) const
G4double GetSandiaCofForMaterial(G4int, G4int) const
static const G4double fZtoAratio[101]
std::vector< G4Material * > G4MaterialTable
static constexpr double cm2
G4double GetPhotoAbsorpCof(G4int i, G4int j) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4double GetWaterEnergyLimit() const
void ComputeMatSandiaMatrixPAI()
std::vector< G4Element * > G4ElementVector
G4OrderedTable * fMatSandiaMatrixPAI
static const G4double fSandiaTable[981][5]
G4GLOB_DLL std::ostream G4cout
static G4int fCumulInterval[101]
static const G4double fH2OlowerI1[23][5]
G4OrderedTable * fMatSandiaMatrix
static constexpr double amu
G4double GetSandiaMatTable(G4int, G4int) const
static const G4double funitc[5]
static constexpr double keV
G4double GetWaterCofForMaterial(G4int, G4int) const
size_t GetNumberOfElements() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double GetDensity() const