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 #include "G4IonisParamMat.hh"
00045 #include "G4Material.hh"
00046 #include "G4DensityEffectData.hh"
00047 #include "G4NistManager.hh"
00048 #include "G4Pow.hh"
00049 #include "G4PhysicalConstants.hh"
00050 #include "G4SystemOfUnits.hh"
00051
00052 G4DensityEffectData* G4IonisParamMat::fDensityData = 0;
00053
00054
00055
00056 G4IonisParamMat::G4IonisParamMat(G4Material* material)
00057 : fMaterial(material)
00058 {
00059 fBirks = 0.;
00060 fMeanEnergyPerIon = 0.0;
00061 twoln10 = 2.*G4Pow::GetInstance()->logZ(10);
00062
00063
00064 fCdensity = 0.0;
00065 fD0density = 0.0;
00066 fAdjustmentFactor = 1.0;
00067 if(!fDensityData) { fDensityData = new G4DensityEffectData(); }
00068
00069
00070 ComputeMeanParameters();
00071 ComputeDensityEffect();
00072 ComputeFluctModel();
00073 ComputeIonParameters();
00074 }
00075
00076
00077
00078
00079
00080
00081 G4IonisParamMat::G4IonisParamMat(__void__&)
00082 : fMaterial(0), fShellCorrectionVector(0)
00083 {
00084 fMeanExcitationEnergy = 0.0;
00085 fLogMeanExcEnergy = 0.0;
00086 fTaul = 0.0;
00087 fCdensity = 0.0;
00088 fMdensity = 0.0;
00089 fAdensity = 0.0;
00090 fX0density = 0.0;
00091 fX1density = 0.0;
00092 fD0density = 0.0;
00093 fPlasmaEnergy = 0.0;
00094 fAdjustmentFactor = 0.0;
00095 fF1fluct = 0.0;
00096 fF2fluct = 0.0;
00097 fEnergy1fluct = 0.0;
00098 fLogEnergy1fluct = 0.0;
00099 fEnergy2fluct = 0.0;
00100 fLogEnergy2fluct = 0.0;
00101 fEnergy0fluct = 0.0;
00102 fRateionexcfluct = 0.0;
00103 fZeff = 0.0;
00104 fFermiEnergy = 0.0;
00105 fLfactor = 0.0;
00106 fInvA23 = 0.0;
00107 fBirks = 0.0;
00108 fMeanEnergyPerIon = 0.0;
00109 twoln10 = 2.*G4Pow::GetInstance()->logZ(10);
00110 }
00111
00112
00113
00114 G4IonisParamMat::~G4IonisParamMat()
00115 {
00116 if (fShellCorrectionVector) { delete [] fShellCorrectionVector; }
00117 if (fDensityData) { delete fDensityData; }
00118 fDensityData = 0;
00119 fShellCorrectionVector = 0;
00120 }
00121
00122
00123
00124 void G4IonisParamMat::ComputeMeanParameters()
00125 {
00126
00127 fTaul = (*(fMaterial->GetElementVector()))[0]->GetIonisation()->GetTaul();
00128
00129 fMeanExcitationEnergy = 0.;
00130 fLogMeanExcEnergy = 0.;
00131
00132 size_t nElements = fMaterial->GetNumberOfElements();
00133 const G4ElementVector* elmVector = fMaterial->GetElementVector();
00134 const G4double* nAtomsPerVolume = fMaterial->GetVecNbOfAtomsPerVolume();
00135
00136 const G4String ch = fMaterial->GetChemicalFormula();
00137
00138 if(ch != "") { fMeanExcitationEnergy = FindMeanExcitationEnergy(ch); }
00139
00140
00141 if(fMeanExcitationEnergy > 0.0) {
00142 fLogMeanExcEnergy = std::log(fMeanExcitationEnergy);
00143
00144
00145 } else {
00146 for (size_t i=0; i < nElements; i++) {
00147 const G4Element* elm = (*elmVector)[i];
00148 fLogMeanExcEnergy += nAtomsPerVolume[i]*elm->GetZ()
00149 *std::log(elm->GetIonisation()->GetMeanExcitationEnergy());
00150 }
00151 fLogMeanExcEnergy /= fMaterial->GetTotNbOfElectPerVolume();
00152 fMeanExcitationEnergy = std::exp(fLogMeanExcEnergy);
00153 }
00154
00155 fShellCorrectionVector = new G4double[3];
00156
00157 for (G4int j=0; j<=2; j++)
00158 {
00159 fShellCorrectionVector[j] = 0.;
00160
00161 for (size_t k=0; k<nElements; k++) {
00162 fShellCorrectionVector[j] += nAtomsPerVolume[k]
00163 *(((*elmVector)[k])->GetIonisation()->GetShellCorrectionVector())[j];
00164 }
00165 fShellCorrectionVector[j] *= 2.0/fMaterial->GetTotNbOfElectPerVolume();
00166 }
00167 }
00168
00169
00170
00171 G4DensityEffectData* G4IonisParamMat::GetDensityEffectData()
00172 {
00173 return fDensityData;
00174 }
00175
00176
00177
00178 void G4IonisParamMat::ComputeDensityEffect()
00179 {
00180 G4State State = fMaterial->GetState();
00181
00182
00183
00184 G4int idx = fDensityData->GetIndex(fMaterial->GetName());
00185 G4int nelm= fMaterial->GetNumberOfElements();
00186 G4int Z0 = G4int((*(fMaterial->GetElementVector()))[0]->GetZ()+0.5);
00187 if(idx < 0 && 1 == nelm) {
00188 idx = fDensityData->GetElementIndex(Z0, fMaterial->GetState());
00189 }
00190
00191
00192
00193 if(idx >= 0) {
00194
00195
00196
00197
00198
00199
00200 fCdensity = fDensityData->GetCdensity(idx);
00201 fMdensity = fDensityData->GetMdensity(idx);
00202 fAdensity = fDensityData->GetAdensity(idx);
00203 fX0density = fDensityData->GetX0density(idx);
00204 fX1density = fDensityData->GetX1density(idx);
00205 fD0density = fDensityData->GetDelta0density(idx);
00206 fPlasmaEnergy = fDensityData->GetPlasmaEnergy(idx);
00207 fAdjustmentFactor = fDensityData->GetAdjustmentFactor(idx);
00208
00209
00210 const G4Material* bmat = fMaterial->GetBaseMaterial();
00211 if(bmat) {
00212 G4double corr = std::log(bmat->GetDensity()/fMaterial->GetDensity());
00213 fCdensity += corr;
00214 fX0density += corr/twoln10;
00215 fX1density += corr/twoln10;
00216 }
00217
00218 } else {
00219
00220 const G4double Cd2 = 4*pi*hbarc_squared*classic_electr_radius;
00221 fPlasmaEnergy = std::sqrt(Cd2*fMaterial->GetTotNbOfElectPerVolume());
00222
00223
00224
00225 G4int icase;
00226
00227 fCdensity = 1. + 2*std::log(fMeanExcitationEnergy/fPlasmaEnergy);
00228
00229
00230
00231 if ((State == kStateSolid)||(State == kStateLiquid)) {
00232
00233 const G4double E100eV = 100.*eV;
00234 const G4double ClimiS[] = {3.681 , 5.215 };
00235 const G4double X0valS[] = {1.0 , 1.5 };
00236 const G4double X1valS[] = {2.0 , 3.0 };
00237
00238 if(fMeanExcitationEnergy < E100eV) { icase = 0; }
00239 else { icase = 1; }
00240
00241 if(fCdensity < ClimiS[icase]) { fX0density = 0.2; }
00242 else { fX0density = 0.326*fCdensity - X0valS[icase]; }
00243
00244 fX1density = X1valS[icase]; fMdensity = 3.0;
00245
00246
00247 if (1 == nelm && 1 == Z0) {
00248 fX0density = 0.425; fX1density = 2.0; fMdensity = 5.949;
00249 }
00250 }
00251
00252
00253
00254 if (State == kStateGas) {
00255
00256 const G4double ClimiG[] = { 10. , 10.5 , 11. , 11.5 , 12.25 , 13.804};
00257 const G4double X0valG[] = { 1.6 , 1.7 , 1.8 , 1.9 , 2.0 , 2.0 };
00258 const G4double X1valG[] = { 4.0 , 4.0 , 4.0 , 4.0 , 4.0 , 5.0 };
00259
00260 icase = 5;
00261 fX0density = 0.326*fCdensity-2.5 ; fX1density = 5.0 ; fMdensity = 3. ;
00262 while((icase > 0)&&(fCdensity < ClimiG[icase])) { icase-- ; }
00263 fX0density = X0valG[icase]; fX1density = X1valG[icase];
00264
00265
00266 if (1 == nelm && 1 == Z0) {
00267 fX0density = 1.837; fX1density = 3.0; fMdensity = 4.754;
00268 }
00269
00270
00271 if (1 == nelm && 2 == Z0) {
00272 fX0density = 2.191; fX1density = 3.0; fMdensity = 3.297;
00273 }
00274 }
00275 }
00276
00277
00278
00279
00280
00281
00282 if (State == kStateGas) {
00283 G4double Density = fMaterial->GetDensity();
00284 G4double Pressure = fMaterial->GetPressure();
00285 G4double Temp = fMaterial->GetTemperature();
00286
00287 G4double DensitySTP = Density*STP_Pressure*Temp/(Pressure*STP_Temperature);
00288
00289 G4double ParCorr = std::log(Density/DensitySTP);
00290
00291 fCdensity -= ParCorr;
00292 fX0density -= ParCorr/twoln10;
00293 fX1density -= ParCorr/twoln10;
00294 }
00295
00296
00297 if(0.0 == fD0density) {
00298 G4double Xa = fCdensity/twoln10;
00299 fAdensity = twoln10*(Xa-fX0density)
00300 /std::pow((fX1density-fX0density),fMdensity);
00301 }
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314 }
00315
00316
00317
00318 void G4IonisParamMat::ComputeFluctModel()
00319 {
00320
00321
00322 G4double Zeff = 0.;
00323 for (size_t i=0;i<fMaterial->GetNumberOfElements();i++) {
00324 Zeff += (fMaterial->GetFractionVector())[i]
00325 *((*(fMaterial->GetElementVector()))[i]->GetZ());
00326 }
00327 if (Zeff > 2.) { fF2fluct = 2./Zeff; }
00328 else { fF2fluct = 0.; }
00329
00330 fF1fluct = 1. - fF2fluct;
00331 fEnergy2fluct = 10.*Zeff*Zeff*eV;
00332 fLogEnergy2fluct = std::log(fEnergy2fluct);
00333 fLogEnergy1fluct = (fLogMeanExcEnergy - fF2fluct*fLogEnergy2fluct)
00334 /fF1fluct;
00335 fEnergy1fluct = std::exp(fLogEnergy1fluct);
00336 fEnergy0fluct = 10.*eV;
00337 fRateionexcfluct = 0.4;
00338 }
00339
00340
00341
00342 void G4IonisParamMat::ComputeIonParameters()
00343 {
00344
00345 const G4ElementVector* theElementVector = fMaterial->GetElementVector() ;
00346 const G4double* theAtomicNumDensityVector =
00347 fMaterial->GetAtomicNumDensityVector() ;
00348 const G4int NumberOfElements = fMaterial->GetNumberOfElements() ;
00349
00350
00351
00352 G4double z(0.0), vF(0.0), lF(0.0), norm(0.0), a23(0.0);
00353
00354 if( 1 == NumberOfElements ) {
00355 const G4Element* element = (*theElementVector)[0];
00356 z = element->GetZ();
00357 vF= element->GetIonisation()->GetFermiVelocity();
00358 lF= element->GetIonisation()->GetLFactor();
00359 a23 = 1.0/G4Pow::GetInstance()->Z23(G4int(element->GetN()));
00360
00361 } else {
00362 for (G4int iel=0; iel<NumberOfElements; iel++)
00363 {
00364 const G4Element* element = (*theElementVector)[iel];
00365 const G4double weight = theAtomicNumDensityVector[iel];
00366 norm += weight ;
00367 z += element->GetZ() * weight;
00368 vF += element->GetIonisation()->GetFermiVelocity() * weight;
00369 lF += element->GetIonisation()->GetLFactor() * weight;
00370 a23 += weight/G4Pow::GetInstance()->Z23(G4int(element->GetN()));
00371 }
00372 z /= norm;
00373 vF /= norm;
00374 lF /= norm;
00375 a23 /= norm;
00376 }
00377 fZeff = z;
00378 fLfactor = lF;
00379 fFermiEnergy = 25.*keV*vF*vF;
00380 fInvA23 = a23;
00381 }
00382
00383
00384
00385 void G4IonisParamMat::SetMeanExcitationEnergy(G4double value)
00386 {
00387 if(value == fMeanExcitationEnergy || value <= 0.0) { return; }
00388 if (G4NistManager::Instance()->GetVerbose() > 1) {
00389 G4cout << "G4Material: Mean excitation energy is changed for "
00390 << fMaterial->GetName()
00391 << " Iold= " << fMeanExcitationEnergy/eV
00392 << "eV; Inew= " << value/eV << " eV;"
00393 << G4endl;
00394 }
00395
00396 fMeanExcitationEnergy = value;
00397
00398
00399 G4double newlog = std::log(value);
00400 G4double corr = 2*(newlog - fLogMeanExcEnergy);
00401 fCdensity += corr;
00402 fX0density += corr/twoln10;
00403 fX1density += corr/twoln10;
00404
00405
00406 fLogMeanExcEnergy = newlog;
00407 ComputeFluctModel();
00408 }
00409
00410
00411
00412 G4double G4IonisParamMat::FindMeanExcitationEnergy(const G4String& chFormula)
00413 {
00414
00415
00416
00417
00418 const size_t numberOfMolecula = 54;
00419 static G4String name[numberOfMolecula] = {
00420
00421 "NH_3", "C_4H_10", "CO_2", "C_2H_6", "C_7H_16",
00422 "C_6H_14", "CH_4", "NO", "N_2O", "C_8H_18",
00423 "C_5H_12", "C_3H_8", "H_2O-Gas",
00424
00425
00426 "C_3H_6O", "C_6H_5NH_2", "C_6H_6", "C_4H_9OH", "CCl_4",
00427 "C_6H_5Cl", "CHCl_3", "C_6H_12", "C_6H_4Cl_2", "C_4Cl_2H_8O",
00428 "C_2Cl_2H_4", "(C_2H_5)_2O", "C_2H_5OH", "C_3H_5(OH)_3","C_7H_16",
00429 "C_6H_14", "CH_3OH", "C_6H_5NO_2","C_5H_12", "C_3H_7OH",
00430 "C_5H_5N", "C_8H_8", "C_2Cl_4", "C_7H_8", "C_2Cl_3H",
00431 "H_2O", "C_8H_10",
00432
00433
00434 "C_5H_5N_5", "C_5H_5N_5O", "(C_6H_11NO)-nylon", "C_25H_52",
00435 "(C_2H_4)-Polyethylene", "(C_5H_8O-2)-Polymethil_Methacrylate",
00436 "(C_8H_8)-Polystyrene", "A-150-tissue", "Al_2O_3", "CaF_2",
00437 "LiF", "Photo_Emulsion", "(C_2F_4)-Teflon", "SiO_2"
00438 };
00439
00440 static G4double meanExcitation[numberOfMolecula] = {
00441
00442 53.7, 48.3, 85.0, 45.4, 49.2,
00443 49.1, 41.7, 87.8, 84.9, 49.5,
00444 48.2, 47.1, 71.6,
00445
00446 64.2, 66.2, 63.4, 59.9, 166.3,
00447 89.1, 156.0, 56.4, 106.5, 103.3,
00448 111.9, 60.0, 62.9, 72.6, 54.4,
00449 54.0, 67.6, 75.8, 53.6, 61.1,
00450 66.2, 64.0, 159.2, 62.5, 148.1,
00451 75.0, 61.8,
00452
00453 71.4, 75.0, 63.9, 48.3, 57.4,
00454 74.0, 68.7, 65.1, 145.2, 166.,
00455 94.0, 331.0, 99.1, 139.2
00456 };
00457
00458 G4double x = fMeanExcitationEnergy;
00459
00460 for(size_t i=0; i<numberOfMolecula; i++) {
00461 if(chFormula == name[i]) {
00462 x = meanExcitation[i]*eV;
00463 break;
00464 }
00465 }
00466 return x;
00467 }
00468
00469
00470
00471 G4IonisParamMat::G4IonisParamMat(const G4IonisParamMat& right)
00472 {
00473 fShellCorrectionVector = 0;
00474 fMaterial = 0;
00475 *this = right;
00476 }
00477
00478
00479
00480 G4IonisParamMat& G4IonisParamMat::operator=(const G4IonisParamMat& right)
00481 {
00482 if (this != &right)
00483 {
00484 fMaterial = right.fMaterial;
00485 fMeanExcitationEnergy = right.fMeanExcitationEnergy;
00486 fLogMeanExcEnergy = right.fLogMeanExcEnergy;
00487 if(fShellCorrectionVector){ delete [] fShellCorrectionVector; }
00488 fShellCorrectionVector = new G4double[3];
00489 fShellCorrectionVector[0] = right.fShellCorrectionVector[0];
00490 fShellCorrectionVector[1] = right.fShellCorrectionVector[1];
00491 fShellCorrectionVector[2] = right.fShellCorrectionVector[2];
00492 fTaul = right.fTaul;
00493 fCdensity = right.fCdensity;
00494 fMdensity = right.fMdensity;
00495 fAdensity = right.fAdensity;
00496 fX0density = right.fX0density;
00497 fX1density = right.fX1density;
00498 fD0density = right.fD0density;
00499 fPlasmaEnergy = right.fPlasmaEnergy;
00500 fAdjustmentFactor = right.fAdjustmentFactor;
00501 fF1fluct = right.fF1fluct;
00502 fF2fluct = right.fF2fluct;
00503 fEnergy1fluct = right.fEnergy1fluct;
00504 fLogEnergy1fluct = right.fLogEnergy1fluct;
00505 fEnergy2fluct = right.fEnergy2fluct;
00506 fLogEnergy2fluct = right.fLogEnergy2fluct;
00507 fEnergy0fluct = right.fEnergy0fluct;
00508 fRateionexcfluct = right.fRateionexcfluct;
00509 fZeff = right.fZeff;
00510 fFermiEnergy = right.fFermiEnergy;
00511 fLfactor = right.fLfactor;
00512 fInvA23 = right.fInvA23;
00513 fBirks = right.fBirks;
00514 fMeanEnergyPerIon = right.fMeanEnergyPerIon;
00515 fDensityData = right.fDensityData;
00516 twoln10 = right.twoln10;
00517 }
00518 return *this;
00519 }
00520
00521
00522
00523 G4int G4IonisParamMat::operator==(const G4IonisParamMat& right) const
00524 {
00525 return (this == (G4IonisParamMat*) &right);
00526 }
00527
00528
00529
00530 G4int G4IonisParamMat::operator!=(const G4IonisParamMat& right) const
00531 {
00532 return (this != (G4IonisParamMat*) &right);
00533 }
00534
00535
00536