00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046 #include "G4SandiaTable.hh"
00047 #include "G4StaticSandiaData.hh"
00048 #include "G4Material.hh"
00049 #include "G4MaterialTable.hh"
00050 #include "G4PhysicalConstants.hh"
00051 #include "G4SystemOfUnits.hh"
00052
00053 G4int G4SandiaTable::fCumulInterval[101] = {0};
00054 G4double G4SandiaTable::fSandiaCofPerAtom[4] = {0.0};
00055 G4double const G4SandiaTable::funitc[5] = {keV,
00056 cm2*keV/g,
00057 cm2*keV*keV/g,
00058 cm2*keV*keV*keV/g,
00059 cm2*keV*keV*keV*keV/g};
00060
00061
00062
00063 G4SandiaTable::G4SandiaTable(G4Material* material)
00064 : fMaterial(material)
00065 {
00066 fMatSandiaMatrix = 0;
00067 fMatSandiaMatrixPAI = 0;
00068 fPhotoAbsorptionCof = 0;
00069
00070 fMatNbOfIntervals = 0;
00071
00072
00073 fMaxInterval = 0;
00074 fVerbose = 0;
00075
00076
00077
00078 fCumulInterval[0] = 1;
00079
00080 for (G4int Z=1; Z<101; ++Z) {
00081 fCumulInterval[Z] = fCumulInterval[Z-1] + fNbOfIntervals[Z];
00082 }
00083
00084
00085 fnulcof[0] = fnulcof[1] = fnulcof[2] = fnulcof[3] = 0.;
00086
00087 fMaxInterval = 0;
00088
00089
00090 ComputeMatSandiaMatrix();
00091
00092
00093 }
00094
00095
00096
00097
00098
00099
00100 G4SandiaTable::G4SandiaTable(__void__&)
00101 : fMaterial(0),fMatSandiaMatrix(0),fMatSandiaMatrixPAI(0),fPhotoAbsorptionCof(0)
00102 {
00103 fnulcof[0] = fnulcof[1] = fnulcof[2] = fnulcof[3] = 0.;
00104 fMaxInterval = 0;
00105 fMatNbOfIntervals = 0;
00106 fVerbose = 0;
00107 }
00108
00109
00110
00111 G4SandiaTable::~G4SandiaTable()
00112 {
00113 if(fMatSandiaMatrix) {
00114 fMatSandiaMatrix->clearAndDestroy();
00115 delete fMatSandiaMatrix;
00116 }
00117 if(fMatSandiaMatrixPAI) {
00118 fMatSandiaMatrixPAI->clearAndDestroy();
00119 delete fMatSandiaMatrixPAI;
00120 }
00121 if(fPhotoAbsorptionCof)
00122 {
00123 delete [] fPhotoAbsorptionCof;
00124 }
00125 }
00126
00127
00128
00129 G4double*
00130 G4SandiaTable::GetSandiaCofPerAtom(G4int Z, G4double energy)
00131 {
00132 G4double Emin = fSandiaTable[fCumulInterval[Z-1]][0]*keV;
00133 G4double Iopot = fIonizationPotentials[Z]*eV;
00134 if (Iopot > Emin) Emin = Iopot;
00135
00136 G4int interval = fNbOfIntervals[Z] - 1;
00137 G4int row = fCumulInterval[Z-1] + interval;
00138 while ((interval>0) && (energy<fSandiaTable[row][0]*keV)) {
00139 --interval;
00140 row = fCumulInterval[Z-1] + interval;
00141 }
00142 if (energy >= Emin)
00143 {
00144 G4double AoverAvo = Z*amu/fZtoAratio[Z];
00145
00146 fSandiaCofPerAtom[0]=AoverAvo*funitc[1]*fSandiaTable[row][1];
00147 fSandiaCofPerAtom[1]=AoverAvo*funitc[2]*fSandiaTable[row][2];
00148 fSandiaCofPerAtom[2]=AoverAvo*funitc[3]*fSandiaTable[row][3];
00149 fSandiaCofPerAtom[3]=AoverAvo*funitc[4]*fSandiaTable[row][4];
00150 }
00151 else
00152 {
00153 fSandiaCofPerAtom[0] = fSandiaCofPerAtom[1] = fSandiaCofPerAtom[2] =
00154 fSandiaCofPerAtom[3] = 0.;
00155 }
00156 return fSandiaCofPerAtom;
00157 }
00158
00159
00160
00161 G4double G4SandiaTable::GetZtoA(G4int Z)
00162 {
00163 assert (Z>0 && Z<101);
00164 return fZtoAratio[Z];
00165 }
00166
00167
00168
00169 void G4SandiaTable::ComputeMatSandiaMatrix()
00170 {
00171
00172 const G4int NbElm = fMaterial->GetNumberOfElements();
00173 const G4ElementVector* ElementVector = fMaterial->GetElementVector();
00174
00175 G4int* Z = new G4int[NbElm];
00176
00177
00178
00179 G4int MaxIntervals = 0;
00180 G4int elm;
00181
00182 for ( elm = 0; elm < NbElm; elm++ )
00183 {
00184 Z[elm] = (G4int)(*ElementVector)[elm]->GetZ();
00185 MaxIntervals += fNbOfIntervals[Z[elm]];
00186 }
00187
00188
00189
00190
00191 G4double* tmp1 = new G4double[MaxIntervals];
00192 G4double IonizationPot;
00193 G4int interval1 = 0;
00194
00195 for ( elm = 0; elm < NbElm; elm++ )
00196 {
00197 IonizationPot = GetIonizationPot(Z[elm]);
00198
00199 for (G4int row = fCumulInterval[ Z[elm]-1 ]; row < fCumulInterval[Z[elm]]; row++ )
00200 {
00201 tmp1[interval1++] = std::max(fSandiaTable[row][0]*keV,IonizationPot);
00202 }
00203 }
00204
00205
00206
00207
00208 G4double* tmp2 = new G4double[MaxIntervals];
00209 G4double Emin;
00210 G4int interval2 = 0;
00211
00212 do
00213 {
00214 Emin = DBL_MAX;
00215
00216 for ( G4int i1 = 0; i1 < MaxIntervals; i1++ )
00217 {
00218 if (tmp1[i1] < Emin) Emin = tmp1[i1];
00219 }
00220 if (Emin < DBL_MAX) tmp2[interval2++] = Emin;
00221
00222 for ( G4int j1 = 0; j1 < MaxIntervals; j1++ )
00223 {
00224 if (tmp1[j1] <= Emin) tmp1[j1] = DBL_MAX;
00225 }
00226 } while (Emin < DBL_MAX);
00227
00228
00229
00230 fMatSandiaMatrix = new G4OrderedTable();
00231 G4int interval;
00232
00233 for (interval = 0; interval < interval2; interval++ )
00234 {
00235 fMatSandiaMatrix->push_back( new G4DataVector(5,0.) );
00236 }
00237
00238
00239
00240 const G4double* NbOfAtomsPerVolume = fMaterial->GetVecNbOfAtomsPerVolume();
00241
00242 const G4double prec = 1.e-03*eV;
00243 G4double coef, oldsum(0.), newsum(0.);
00244 fMatNbOfIntervals = 0;
00245
00246 for ( interval = 0; interval < interval2; interval++ )
00247 {
00248 Emin = (*(*fMatSandiaMatrix)[fMatNbOfIntervals])[0] = tmp2[interval];
00249
00250 for ( G4int k = 1; k < 5; k++ ) {
00251 (*(*fMatSandiaMatrix)[fMatNbOfIntervals])[k] = 0.;
00252 }
00253 newsum = 0.;
00254
00255 for ( elm = 0; elm < NbElm; elm++ )
00256 {
00257 GetSandiaCofPerAtom(Z[elm], Emin+prec);
00258
00259 for ( G4int j = 1; j < 5; j++ )
00260 {
00261 coef = NbOfAtomsPerVolume[elm]*fSandiaCofPerAtom[j-1];
00262 (*(*fMatSandiaMatrix)[fMatNbOfIntervals])[j] += coef;
00263 newsum += std::fabs(coef);
00264 }
00265 }
00266
00267
00268 if (newsum != oldsum) { oldsum = newsum; fMatNbOfIntervals++;}
00269 }
00270 delete [] Z;
00271 delete [] tmp1;
00272 delete [] tmp2;
00273
00274 if ( fVerbose > 0 && fMaterial->GetName() == "G4_Ar" )
00275 {
00276 G4cout<<"mma, G4SandiaTable::ComputeMatSandiaMatrix(), mat = "
00277 <<fMaterial->GetName()<<G4endl;
00278
00279 for( G4int i = 0; i < fMatNbOfIntervals; i++)
00280 {
00281 G4cout<<i<<"\t"<<GetSandiaCofForMaterial(i,0)/keV<<" keV \t"
00282 <<this->GetSandiaCofForMaterial(i,1)
00283 <<"\t"<<this->GetSandiaCofForMaterial(i,2)
00284 <<"\t"<<this->GetSandiaCofForMaterial(i,3)
00285 <<"\t"<<this->GetSandiaCofForMaterial(i,4)<<G4endl;
00286 }
00287 }
00288 }
00289
00290
00291
00292 void G4SandiaTable::ComputeMatSandiaMatrixPAI()
00293 {
00294 G4int MaxIntervals = 0;
00295 G4int elm, c, i, j, jj, k, k1, k2, c1, n1;
00296
00297 const G4int noElm = fMaterial->GetNumberOfElements();
00298 const G4ElementVector* ElementVector = fMaterial->GetElementVector();
00299 G4int* Z = new G4int[noElm];
00300
00301 for ( elm = 0; elm < noElm; elm++ )
00302 {
00303 Z[elm] = (G4int)(*ElementVector)[elm]->GetZ();
00304 MaxIntervals += fNbOfIntervals[Z[elm]];
00305 }
00306 fMaxInterval = MaxIntervals + 2;
00307
00308 if ( fVerbose > 0 && fMaterial->GetName() == "G4_Ar" )
00309 {
00310 G4cout<<"fMaxInterval = "<<fMaxInterval<<G4endl;
00311 }
00312 G4double* fPhotoAbsorptionCof0 = new G4double[fMaxInterval];
00313 G4double* fPhotoAbsorptionCof1 = new G4double[fMaxInterval];
00314 G4double* fPhotoAbsorptionCof2 = new G4double[fMaxInterval];
00315 G4double* fPhotoAbsorptionCof3 = new G4double[fMaxInterval];
00316 G4double* fPhotoAbsorptionCof4 = new G4double[fMaxInterval];
00317
00318 for( c = 0; c < fMaxInterval; c++ )
00319 {
00320 fPhotoAbsorptionCof0[c] = 0.;
00321 fPhotoAbsorptionCof1[c] = 0.;
00322 fPhotoAbsorptionCof2[c] = 0.;
00323 fPhotoAbsorptionCof3[c] = 0.;
00324 fPhotoAbsorptionCof4[c] = 0.;
00325 }
00326 c = 1;
00327
00328 for(i = 0; i < noElm; i++)
00329 {
00330 G4double I1 = fIonizationPotentials[Z[i]]*keV;
00331 n1 = 1;
00332
00333 for( j = 1; j < Z[i]; j++ ) n1 += fNbOfIntervals[j];
00334
00335 G4int n2 = n1 + fNbOfIntervals[Z[i]];
00336
00337 for( k1 = n1; k1 < n2; k1++ )
00338 {
00339 if( I1 > fSandiaTable[k1][0] )
00340 {
00341 continue;
00342 }
00343 break;
00344 }
00345 G4int flag = 0;
00346
00347 for( c1 = 1; c1 < c; c1++ )
00348 {
00349 if( fPhotoAbsorptionCof0[c1] == I1 )
00350 {
00351 flag = 1;
00352 break;
00353 }
00354 }
00355 if(flag == 0)
00356 {
00357 fPhotoAbsorptionCof0[c] = I1;
00358 c++;
00359 }
00360 for( k2 = k1; k2 < n2; k2++ )
00361 {
00362 flag = 0;
00363
00364 for( c1 = 1; c1 < c; c1++ )
00365 {
00366 if( fPhotoAbsorptionCof0[c1] == fSandiaTable[k2][0] )
00367 {
00368 flag = 1;
00369 break;
00370 }
00371 }
00372 if(flag == 0)
00373 {
00374 fPhotoAbsorptionCof0[c] = fSandiaTable[k2][0];
00375 c++;
00376 }
00377 }
00378 }
00379
00380
00381 for( i = 1; i < c; i++ )
00382 {
00383 for( j = i + 1; j < c; j++ )
00384 {
00385 if( fPhotoAbsorptionCof0[i] > fPhotoAbsorptionCof0[j] )
00386 {
00387 G4double tmp = fPhotoAbsorptionCof0[i];
00388 fPhotoAbsorptionCof0[i] = fPhotoAbsorptionCof0[j];
00389 fPhotoAbsorptionCof0[j] = tmp;
00390 }
00391 }
00392 if ( fVerbose > 0 && fMaterial->GetName() == "G4_Ar" )
00393 {
00394 G4cout<<i<<"\t energy = "<<fPhotoAbsorptionCof0[i]<<G4endl;
00395 }
00396 }
00397 fMaxInterval = c;
00398
00399 const G4double* fractionW = fMaterial->GetFractionVector();
00400
00401 if ( fVerbose > 0 && fMaterial->GetName() == "G4_Ar" )
00402 {
00403 for( i = 0; i < noElm; i++ ) G4cout<<i<<" = elN, fraction = "<<fractionW[i]<<G4endl;
00404 }
00405
00406 for( i = 0; i < noElm; i++ )
00407 {
00408 n1 = 1;
00409 G4double I1 = fIonizationPotentials[Z[i]]*keV;
00410
00411 for( j = 1; j < Z[i]; j++ ) n1 += fNbOfIntervals[j];
00412
00413 G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1;
00414
00415 for(k = n1; k < n2; k++)
00416 {
00417 G4double B1 = fSandiaTable[k][0];
00418 G4double B2 = fSandiaTable[k+1][0];
00419
00420 for(G4int q = 1; q < fMaxInterval-1; q++)
00421 {
00422 G4double E1 = fPhotoAbsorptionCof0[q];
00423 G4double E2 = fPhotoAbsorptionCof0[q+1];
00424
00425 if ( fVerbose > 0 && fMaterial->GetName() == "G4_Ar" )
00426 {
00427 G4cout<<"k = "<<k<<", q = "<<q<<", B1 = "<<B1<<", B2 = "<<B2<<", E1 = "<<E1<<", E2 = "<<E2<<G4endl;
00428 }
00429 if( B1 > E1 || B2 < E2 || E1 < I1 )
00430 {
00431 if ( fVerbose > 0 && fMaterial->GetName() == "G4_Ar" )
00432 {
00433 G4cout<<"continue for: B1 = "<<B1<<", B2 = "<<B2<<", E1 = "<<E1<<", E2 = "<<E2<<G4endl;
00434 }
00435 continue;
00436 }
00437 fPhotoAbsorptionCof1[q] += fSandiaTable[k][1]*fractionW[i];
00438 fPhotoAbsorptionCof2[q] += fSandiaTable[k][2]*fractionW[i];
00439 fPhotoAbsorptionCof3[q] += fSandiaTable[k][3]*fractionW[i];
00440 fPhotoAbsorptionCof4[q] += fSandiaTable[k][4]*fractionW[i];
00441 }
00442 }
00443
00444
00445 fPhotoAbsorptionCof1[fMaxInterval-1] += fSandiaTable[k][1]*fractionW[i];
00446 fPhotoAbsorptionCof2[fMaxInterval-1] += fSandiaTable[k][2]*fractionW[i];
00447 fPhotoAbsorptionCof3[fMaxInterval-1] += fSandiaTable[k][3]*fractionW[i];
00448 fPhotoAbsorptionCof4[fMaxInterval-1] += fSandiaTable[k][4]*fractionW[i];
00449 }
00450 c = 0;
00451
00452 do
00453 {
00454 c++;
00455
00456 if( fPhotoAbsorptionCof1[c] != 0.0 ||
00457 fPhotoAbsorptionCof2[c] != 0.0 ||
00458 fPhotoAbsorptionCof3[c] != 0.0 ||
00459 fPhotoAbsorptionCof4[c] != 0.0 ) continue;
00460
00461 if ( fVerbose > 0 && fMaterial->GetName() == "G4_Ar" )
00462 {
00463 G4cout<<c<<" = number with zero cofs"<<G4endl;
00464 }
00465 for( jj = 2; jj < fMaxInterval; jj++ )
00466 {
00467 fPhotoAbsorptionCof0[jj-1] = fPhotoAbsorptionCof0[jj];
00468 fPhotoAbsorptionCof1[jj-1] = fPhotoAbsorptionCof1[jj];
00469 fPhotoAbsorptionCof2[jj-1] = fPhotoAbsorptionCof2[jj];
00470 fPhotoAbsorptionCof3[jj-1] = fPhotoAbsorptionCof3[jj];
00471 fPhotoAbsorptionCof4[jj-1] = fPhotoAbsorptionCof4[jj];
00472 }
00473 fMaxInterval--;
00474 c--;
00475 }
00476 while( c < fMaxInterval - 1 );
00477
00478
00479
00480 fMatSandiaMatrixPAI = new G4OrderedTable();
00481 G4double density = fMaterial->GetDensity();
00482
00483 for (i = 0; i < fMaxInterval; i++) fMatSandiaMatrixPAI->push_back(new G4DataVector(5,0.));
00484
00485 for (i = 0; i < fMaxInterval; i++)
00486 {
00487 (*(*fMatSandiaMatrixPAI)[i])[0] = fPhotoAbsorptionCof0[i+1];
00488 (*(*fMatSandiaMatrixPAI)[i])[1] = fPhotoAbsorptionCof1[i+1]*density;
00489 (*(*fMatSandiaMatrixPAI)[i])[2] = fPhotoAbsorptionCof2[i+1]*density;
00490 (*(*fMatSandiaMatrixPAI)[i])[3] = fPhotoAbsorptionCof3[i+1]*density;
00491 (*(*fMatSandiaMatrixPAI)[i])[4] = fPhotoAbsorptionCof4[i+1]*density;
00492 }
00493 if ( fVerbose > 0 && fMaterial->GetName() == "G4_Ar" )
00494 {
00495 G4cout<<"mma, G4SandiaTable::ComputeMatSandiaMatrixPAI(), mat = "<<fMaterial->GetName()<<G4endl;
00496
00497 for( i = 0; i < fMaxInterval; i++)
00498 {
00499 G4cout<<i<<"\t"<<GetSandiaMatTablePAI(i,0)/keV<<" keV \t"<<this->GetSandiaMatTablePAI(i,1)
00500 <<"\t"<<this->GetSandiaMatTablePAI(i,2)<<"\t"<<this->GetSandiaMatTablePAI(i,3)
00501 <<"\t"<<this->GetSandiaMatTablePAI(i,4)<<G4endl;
00502 }
00503 }
00504
00505 delete [] Z;
00506 delete [] fPhotoAbsorptionCof0;
00507 delete [] fPhotoAbsorptionCof1;
00508 delete [] fPhotoAbsorptionCof2;
00509 delete [] fPhotoAbsorptionCof3;
00510 delete [] fPhotoAbsorptionCof4;
00511 return;
00512 }
00513
00514
00515
00518
00519
00520
00521
00522
00523 G4SandiaTable::G4SandiaTable(G4int matIndex)
00524 {
00525 fMaterial = 0;
00526 fMatNbOfIntervals = 0;
00527 fMatSandiaMatrix = 0;
00528 fMatSandiaMatrixPAI = 0;
00529 fPhotoAbsorptionCof = 0;
00530
00531
00532 fMaxInterval = 0;
00533 fVerbose = 0;
00534
00535
00536 fnulcof[0] = fnulcof[1] = fnulcof[2] = fnulcof[3] = 0.;
00537
00538 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
00539 G4int numberOfMat = G4Material::GetNumberOfMaterials();
00540
00541 if ( matIndex >= 0 && matIndex < numberOfMat)
00542 {
00543 fMaterial = (*theMaterialTable)[matIndex];
00544 ComputeMatTable();
00545 }
00546 else
00547 {
00548 G4Exception("G4SandiaTable::G4SandiaTable(G4int matIndex)", "mat401",
00549 FatalException, "wrong matIndex");
00550 }
00551 }
00552
00554
00555
00556
00557
00558 void
00559 G4SandiaTable::SandiaSort(G4double** da ,
00560 G4int sz )
00561 {
00562 for(G4int i = 1;i < sz; i++ )
00563 {
00564 for(G4int j = i + 1;j < sz; j++ )
00565 {
00566 if(da[i][0] > da[j][0]) SandiaSwap(da,i,j);
00567 }
00568 }
00569 }
00570
00572
00573
00574
00575
00576 G4int
00577 G4SandiaTable::SandiaIntervals(G4int Z[],
00578 G4int el )
00579 {
00580 G4int c, i, flag = 0, n1 = 1;
00581 G4int j, c1, k1, k2;
00582 G4double I1;
00583 fMaxInterval = 0;
00584
00585 for( i = 0; i < el; i++ ) fMaxInterval += fNbOfIntervals[ Z[i] ];
00586
00587 fMaxInterval += 2;
00588
00589 if( fVerbose > 0 ) G4cout<<"begin sanInt, fMaxInterval = "<<fMaxInterval<<G4endl;
00590
00591 fPhotoAbsorptionCof = new G4double* [fMaxInterval];
00592
00593 for( i = 0; i < fMaxInterval; i++ ) fPhotoAbsorptionCof[i] = new G4double[5];
00594
00595
00596
00597 for( c = 0; c < fMaxInterval; c++ ) fPhotoAbsorptionCof[c][0] = 0.;
00598
00599 c = 1;
00600
00601 for( i = 0; i < el; i++ )
00602 {
00603 I1 = fIonizationPotentials[ Z[i] ]*keV;
00604 n1 = 1;
00605
00606 for( j = 1; j < Z[i]; j++ ) n1 += fNbOfIntervals[j];
00607
00608 G4int n2 = n1 + fNbOfIntervals[Z[i]];
00609
00610 for( k1 = n1; k1 < n2; k1++ )
00611 {
00612 if( I1 > fSandiaTable[k1][0] )
00613 {
00614 continue;
00615 }
00616 break;
00617 }
00618 flag = 0;
00619
00620 for( c1 = 1; c1 < c; c1++ )
00621 {
00622 if( fPhotoAbsorptionCof[c1][0] == I1 )
00623 {
00624 flag = 1;
00625 break;
00626 }
00627 }
00628 if( flag == 0 )
00629 {
00630 fPhotoAbsorptionCof[c][0] = I1;
00631 c++;
00632 }
00633 for( k2 = k1; k2 < n2; k2++ )
00634 {
00635 flag = 0;
00636
00637 for( c1 = 1; c1 < c; c1++ )
00638 {
00639 if( fPhotoAbsorptionCof[c1][0] == fSandiaTable[k2][0] )
00640 {
00641 flag = 1;
00642 break;
00643 }
00644 }
00645 if( flag == 0 )
00646 {
00647 fPhotoAbsorptionCof[c][0] = fSandiaTable[k2][0];
00648 if( fVerbose > 0 ) G4cout<<"sanInt, c = "<<c<<", E_c = "<<fPhotoAbsorptionCof[c][0]<<G4endl;
00649 c++;
00650 }
00651 }
00652 }
00653
00654 SandiaSort(fPhotoAbsorptionCof,c);
00655 fMaxInterval = c;
00656 if( fVerbose > 0 ) G4cout<<"end SanInt, fMaxInterval = "<<fMaxInterval<<G4endl;
00657 return c;
00658 }
00659
00661
00662
00663
00664
00665 G4int
00666 G4SandiaTable::SandiaMixing( G4int Z[],
00667 const G4double fractionW[],
00668 G4int el,
00669 G4int mi )
00670 {
00671 G4int i, j, n1, k, c=1, jj, kk;
00672 G4double I1, B1, B2, E1, E2;
00673
00674 for( i = 0; i < mi; i++ )
00675 {
00676 for( j = 1; j < 5; j++ ) fPhotoAbsorptionCof[i][j] = 0.;
00677 }
00678 for( i = 0; i < el; i++ )
00679 {
00680 n1 = 1;
00681 I1 = fIonizationPotentials[Z[i]]*keV;
00682
00683 for( j = 1; j < Z[i]; j++ ) n1 += fNbOfIntervals[j];
00684
00685 G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1;
00686
00687 for( k = n1; k < n2; k++ )
00688 {
00689 B1 = fSandiaTable[k][0];
00690 B2 = fSandiaTable[k+1][0];
00691
00692 for( c = 1; c < mi-1; c++ )
00693 {
00694 E1 = fPhotoAbsorptionCof[c][0];
00695 E2 = fPhotoAbsorptionCof[c+1][0];
00696
00697 if( B1 > E1 || B2 < E2 || E1 < I1 ) continue;
00698
00699 for( j = 1; j < 5; j++ )
00700 {
00701 fPhotoAbsorptionCof[c][j] += fSandiaTable[k][j]*fractionW[i];
00702 if( fVerbose > 0 )
00703 {
00704 G4cout<<"c="<<c<<"; j="<<j<<"; fST="<<fSandiaTable[k][j]<<"; frW="<<fractionW[i]<<G4endl;
00705 }
00706 }
00707 }
00708 }
00709 for( j = 1; j < 5; j++ )
00710 {
00711 fPhotoAbsorptionCof[mi-1][j] += fSandiaTable[k][j]*fractionW[i];
00712 if( fVerbose > 0 )
00713 {
00714 G4cout<<"mi-1="<<mi-1<<"; j="<<j<<"; fST="<<fSandiaTable[k][j]<<"; frW="<<fractionW[i]<<G4endl;
00715 }
00716 }
00717 }
00718 c = 0;
00719
00720 do
00721 {
00722 c++;
00723
00724 if( fPhotoAbsorptionCof[c][1] != 0.0 ||
00725 fPhotoAbsorptionCof[c][2] != 0.0 ||
00726 fPhotoAbsorptionCof[c][3] != 0.0 ||
00727 fPhotoAbsorptionCof[c][4] != 0.0 ) continue;
00728
00729 for( jj = 2; jj < mi; jj++ )
00730 {
00731 for( kk = 0; kk < 5; kk++ ) fPhotoAbsorptionCof[jj-1][kk] = fPhotoAbsorptionCof[jj][kk];
00732 }
00733 mi--;
00734 c--;
00735 }
00736 while( c < mi - 1 );
00737
00738 if( fVerbose > 0 ) G4cout<<"end SanMix, mi = "<<mi<<G4endl;
00739
00740 return mi;
00741 }
00742
00743
00744
00745 G4int G4SandiaTable::GetMatNbOfIntervals()
00746 {
00747 return fMatNbOfIntervals;
00748 }
00749
00750
00751
00752 G4int G4SandiaTable::GetNbOfIntervals(G4int Z)
00753 {
00754 assert (Z>0 && Z<101);
00755 return fNbOfIntervals[Z];
00756 }
00757
00758
00759
00760 G4double
00761 G4SandiaTable::GetSandiaCofPerAtom(G4int Z, G4int interval, G4int j)
00762 {
00763 assert (Z>0 && Z<101 && interval>=0 && interval<fNbOfIntervals[Z]
00764 && j>=0 && j<5);
00765
00766 G4int row = fCumulInterval[Z-1] + interval;
00767 G4double x = fSandiaTable[row][0]*CLHEP::keV;
00768 if (j > 0) {
00769 x = Z*CLHEP::amu/fZtoAratio[Z]*fSandiaTable[row][j]*funitc[j];
00770 }
00771 return x;
00772 }
00773
00774
00775
00776 G4double
00777 G4SandiaTable::GetSandiaCofForMaterial(G4int interval, G4int j)
00778 {
00779 assert (interval>=0 && interval<fMatNbOfIntervals && j>=0 && j<5);
00780 return ((*(*fMatSandiaMatrix)[interval])[j]);
00781 }
00782
00783
00784
00785 G4double*
00786 G4SandiaTable::GetSandiaCofForMaterial(G4double energy)
00787 {
00788 G4double* x = fnulcof;
00789 if (energy >= (*(*fMatSandiaMatrix)[0])[0]) {
00790
00791 G4int interval = fMatNbOfIntervals - 1;
00792 while ((interval>0)&&(energy<(*(*fMatSandiaMatrix)[interval])[0]))
00793 {interval--;}
00794 x = &((*(*fMatSandiaMatrix)[interval])[1]);
00795 }
00796 return x;
00797 }
00798
00799
00800
00801 G4double
00802 G4SandiaTable::GetSandiaMatTable(G4int interval, G4int j)
00803 {
00804 assert (interval >= 0 && interval < fMaxInterval && j >= 0 && j < 5 );
00805 return ((*(*fMatSandiaMatrix)[interval])[j])*funitc[j];
00806 }
00807
00808
00809
00810 G4double
00811 G4SandiaTable::GetSandiaCofForMaterialPAI(G4int interval, G4int j)
00812 {
00813 assert (interval>=0 && interval<fMatNbOfIntervals && j>=0 && j<5);
00814 if(!fMatSandiaMatrixPAI) ComputeMatSandiaMatrixPAI();
00815 return ((*(*fMatSandiaMatrixPAI)[interval])[j]);
00816 }
00817
00818
00819
00820 G4double*
00821 G4SandiaTable::GetSandiaCofForMaterialPAI(G4double energy)
00822 {
00823 if(!fMatSandiaMatrixPAI) ComputeMatSandiaMatrixPAI();
00824 G4double* x = fnulcof;
00825 if (energy >= (*(*fMatSandiaMatrixPAI)[0])[0]) {
00826
00827 G4int interval = fMatNbOfIntervals - 1;
00828 while ((interval>0)&&(energy<(*(*fMatSandiaMatrixPAI)[interval])[0]))
00829 {interval--;}
00830 x = &((*(*fMatSandiaMatrixPAI)[interval])[1]);
00831 }
00832 return x;
00833 }
00834
00835
00836
00837 G4double
00838 G4SandiaTable::GetSandiaMatTablePAI(G4int interval, G4int j)
00839 {
00840 assert (interval >= 0 && interval < fMaxInterval && j >= 0 && j < 5 );
00841 if(!fMatSandiaMatrixPAI) { ComputeMatSandiaMatrixPAI(); }
00842 return ((*(*fMatSandiaMatrixPAI)[interval])[j])*funitc[j];
00843 }
00844
00845
00846
00847 G4double
00848 G4SandiaTable::GetIonizationPot(G4int Z)
00849 {
00850 assert (Z>0 && Z<101);
00851 return fIonizationPotentials[Z]*CLHEP::eV;
00852 }
00853
00854
00855
00856 G4OrderedTable*
00857 G4SandiaTable::GetSandiaMatrixPAI()
00858 {
00859 if(!fMatSandiaMatrixPAI) { ComputeMatSandiaMatrixPAI(); }
00860 return fMatSandiaMatrixPAI;
00861 }
00862
00864
00865
00866
00867
00868 void G4SandiaTable::ComputeMatTable()
00869 {
00870 G4int MaxIntervals = 0;
00871 G4int elm, c, i, j, jj, k, kk, k1, k2, c1, n1;
00872
00873 const G4int noElm = fMaterial->GetNumberOfElements();
00874 const G4ElementVector* ElementVector = fMaterial->GetElementVector();
00875 G4int* Z = new G4int[noElm];
00876
00877 for (elm = 0; elm<noElm; elm++)
00878 {
00879 Z[elm] = (G4int)(*ElementVector)[elm]->GetZ();
00880 MaxIntervals += fNbOfIntervals[Z[elm]];
00881 }
00882 fMaxInterval = 0;
00883
00884 for(i = 0; i < noElm; i++) fMaxInterval += fNbOfIntervals[Z[i]];
00885
00886 fMaxInterval += 2;
00887
00888
00889
00890 fPhotoAbsorptionCof = new G4double* [fMaxInterval];
00891
00892 for(i = 0; i < fMaxInterval; i++)
00893 {
00894 fPhotoAbsorptionCof[i] = new G4double[5];
00895 }
00896
00897
00898
00899 for(c = 0; c < fMaxInterval; c++)
00900 {
00901 fPhotoAbsorptionCof[c][0] = 0.;
00902 }
00903 c = 1;
00904
00905 for(i = 0; i < noElm; i++)
00906 {
00907 G4double I1 = fIonizationPotentials[Z[i]]*keV;
00908 n1 = 1;
00909
00910 for(j = 1; j < Z[i]; j++)
00911 {
00912 n1 += fNbOfIntervals[j];
00913 }
00914 G4int n2 = n1 + fNbOfIntervals[Z[i]];
00915
00916 for(k1 = n1; k1 < n2; k1++)
00917 {
00918 if(I1 > fSandiaTable[k1][0])
00919 {
00920 continue;
00921 }
00922 break;
00923 }
00924 G4int flag = 0;
00925
00926 for(c1 = 1; c1 < c; c1++)
00927 {
00928 if(fPhotoAbsorptionCof[c1][0] == I1)
00929 {
00930 flag = 1;
00931 break;
00932 }
00933 }
00934 if(flag == 0)
00935 {
00936 fPhotoAbsorptionCof[c][0] = I1;
00937 c++;
00938 }
00939 for(k2 = k1; k2 < n2; k2++)
00940 {
00941 flag = 0;
00942
00943 for(c1 = 1; c1 < c; c1++)
00944 {
00945 if(fPhotoAbsorptionCof[c1][0] == fSandiaTable[k2][0])
00946 {
00947 flag = 1;
00948 break;
00949 }
00950 }
00951 if(flag == 0)
00952 {
00953 fPhotoAbsorptionCof[c][0] = fSandiaTable[k2][0];
00954 c++;
00955 }
00956 }
00957 }
00958
00959 SandiaSort(fPhotoAbsorptionCof,c);
00960 fMaxInterval = c;
00961
00962 const G4double* fractionW = fMaterial->GetFractionVector();
00963
00964 for(i = 0; i < fMaxInterval; i++)
00965 {
00966 for(j = 1; j < 5; j++) fPhotoAbsorptionCof[i][j] = 0.;
00967 }
00968 for(i = 0; i < noElm; i++)
00969 {
00970 n1 = 1;
00971 G4double I1 = fIonizationPotentials[Z[i]]*keV;
00972
00973 for(j = 1; j < Z[i]; j++)
00974 {
00975 n1 += fNbOfIntervals[j];
00976 }
00977 G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1;
00978
00979 for(k = n1; k < n2; k++)
00980 {
00981 G4double B1 = fSandiaTable[k][0];
00982 G4double B2 = fSandiaTable[k+1][0];
00983 for(G4int q = 1; q < fMaxInterval-1; q++)
00984 {
00985 G4double E1 = fPhotoAbsorptionCof[q][0];
00986 G4double E2 = fPhotoAbsorptionCof[q+1][0];
00987 if(B1 > E1 || B2 < E2 || E1 < I1)
00988 {
00989 continue;
00990 }
00991 for(j = 1; j < 5; j++)
00992 {
00993 fPhotoAbsorptionCof[q][j] += fSandiaTable[k][j]*fractionW[i];
00994 }
00995 }
00996 }
00997 for(j = 1; j < 5; j++)
00998 {
00999 fPhotoAbsorptionCof[fMaxInterval-1][j] += fSandiaTable[k][j]*fractionW[i];
01000 }
01001 }
01002
01003 c = 0;
01004
01005 do
01006 {
01007 c++;
01008
01009 if( fPhotoAbsorptionCof[c][1] != 0.0 ||
01010 fPhotoAbsorptionCof[c][2] != 0.0 ||
01011 fPhotoAbsorptionCof[c][3] != 0.0 ||
01012 fPhotoAbsorptionCof[c][4] != 0.0 ) continue;
01013
01014 for(jj = 2; jj < fMaxInterval; jj++)
01015 {
01016 for(kk = 0; kk < 5; kk++)
01017 {
01018 fPhotoAbsorptionCof[jj-1][kk]= fPhotoAbsorptionCof[jj][kk];
01019 }
01020 }
01021 fMaxInterval--;
01022 c--;
01023 }
01024 while(c < fMaxInterval - 1);
01025
01026
01027
01028 fMaxInterval--;
01029
01030 fMatSandiaMatrix = new G4OrderedTable();
01031
01032 for (i = 0; i < fMaxInterval; i++)
01033 {
01034 fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
01035 }
01036 for ( i = 0; i < fMaxInterval; i++ )
01037 {
01038 for( j = 0; j < 5; j++ )
01039 {
01040 (*(*fMatSandiaMatrix)[i])[j] = fPhotoAbsorptionCof[i+1][j];
01041 }
01042 }
01043 fMatNbOfIntervals = fMaxInterval;
01044
01045 if ( fVerbose > 0 )
01046 {
01047 G4cout<<"vmg, G4SandiaTable::ComputeMatTable(), mat = "<<fMaterial->GetName()<<G4endl;
01048
01049 for ( i = 0; i < fMaxInterval; i++ )
01050 {
01051
01052
01053
01054
01055 G4cout<<i<<"\t"<<GetSandiaCofForMaterial(i,0)/keV<<" keV \t"<<this->GetSandiaCofForMaterial(i,1)
01056 <<"\t"<<this->GetSandiaCofForMaterial(i,2)<<"\t"<<this->GetSandiaCofForMaterial(i,3)
01057 <<"\t"<<this->GetSandiaCofForMaterial(i,4)<<G4endl;
01058 }
01059 }
01060 delete [] Z;
01061 return;
01062 }
01063
01064
01065
01067
01068