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 #include "G4LivermorePhotoElectricModel.hh"
00037 #include "G4SystemOfUnits.hh"
00038 #include "G4LossTableManager.hh"
00039 #include "G4Electron.hh"
00040 #include "G4Gamma.hh"
00041 #include "G4ParticleChangeForGamma.hh"
00042 #include "G4CrossSectionHandler.hh"
00043 #include "G4LPhysicsFreeVector.hh"
00044 #include "G4VAtomDeexcitation.hh"
00045 #include "G4SauterGavrilaAngularDistribution.hh"
00046 #include "G4AtomicShell.hh"
00047
00048
00049
00050 using namespace std;
00051
00052
00053
00054 G4LivermorePhotoElectricModel::G4LivermorePhotoElectricModel(
00055 const G4String& nam)
00056 : G4VEmModel(nam),fParticleChange(0),maxZ(99),
00057 nShellLimit(100),fDeexcitationActive(false),isInitialised(false),
00058 fAtomDeexcitation(0)
00059 {
00060 verboseLevel= 0;
00061
00062
00063
00064
00065
00066
00067
00068 theGamma = G4Gamma::Gamma();
00069 theElectron = G4Electron::Electron();
00070
00071
00072 SetAngularDistribution(new G4SauterGavrilaAngularDistribution());
00073
00074 for(G4int i=0; i<maxZ; ++i) {
00075 fCrossSection[i] = 0;
00076 fCrossSectionLE[i] = 0;
00077 fNShells[i] = 0;
00078 fNShellsUsed[i] = 0;
00079 }
00080
00081 if(verboseLevel>0) {
00082 G4cout << "Livermore PhotoElectric is constructed "
00083 << " nShellLimit= " << nShellLimit << G4endl;
00084 }
00085
00086
00087 SetDeexcitationFlag(true);
00088 }
00089
00090
00091
00092 G4LivermorePhotoElectricModel::~G4LivermorePhotoElectricModel()
00093 {
00094 for(G4int i=0; i<maxZ; ++i) {
00095 delete fCrossSection[i];
00096 delete fCrossSectionLE[i];
00097 }
00098 }
00099
00100
00101
00102 void
00103 G4LivermorePhotoElectricModel::Initialise(const G4ParticleDefinition*,
00104 const G4DataVector&)
00105 {
00106 if (verboseLevel > 2) {
00107 G4cout << "Calling G4LivermorePhotoElectricModel::Initialise()" << G4endl;
00108 }
00109
00110 char* path = getenv("G4LEDATA");
00111
00112 G4ProductionCutsTable* theCoupleTable =
00113 G4ProductionCutsTable::GetProductionCutsTable();
00114 G4int numOfCouples = theCoupleTable->GetTableSize();
00115
00116 for(G4int i=0; i<numOfCouples; ++i)
00117 {
00118 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
00119 const G4Material* material = couple->GetMaterial();
00120 const G4ElementVector* theElementVector = material->GetElementVector();
00121 G4int nelm = material->GetNumberOfElements();
00122
00123 for (G4int j=0; j<nelm; ++j)
00124 {
00125 G4int Z = (G4int)(*theElementVector)[j]->GetZ();
00126 if(Z < 1) { Z = 1; }
00127 else if(Z > maxZ) { Z = maxZ; }
00128 if(!fCrossSection[Z]) { ReadData(Z, path); }
00129 }
00130 }
00131
00132 if (verboseLevel > 2) {
00133 G4cout << "Loaded cross section files for LivermorePhotoElectric model"
00134 << G4endl;
00135 }
00136 if(!isInitialised) {
00137 isInitialised = true;
00138 fParticleChange = GetParticleChangeForGamma();
00139
00140 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
00141 }
00142 fDeexcitationActive = false;
00143 if(fAtomDeexcitation) {
00144 fDeexcitationActive = fAtomDeexcitation->IsFluoActive();
00145 }
00146
00147 if (verboseLevel > 0) {
00148 G4cout << "LivermorePhotoElectric model is initialized " << G4endl
00149 << G4endl;
00150 }
00151 }
00152
00153
00154
00155 G4double G4LivermorePhotoElectricModel::ComputeCrossSectionPerAtom(
00156 const G4ParticleDefinition*,
00157 G4double energy,
00158 G4double ZZ, G4double,
00159 G4double, G4double)
00160 {
00161 if (verboseLevel > 3) {
00162 G4cout << "G4LivermorePhotoElectricModel::Calling ComputeCrossSectionPerAtom()"
00163 << " Z= " << ZZ << " R(keV)= " << energy/keV << G4endl;
00164 }
00165 G4double cs = 0.0;
00166 G4double gammaEnergy = energy;
00167
00168 G4int Z = G4lrint(ZZ);
00169 if(Z < 1 || Z >= maxZ) { return cs; }
00170
00171
00172 if(!fCrossSection[Z]) {
00173 char* path = getenv("G4LEDATA");
00174 ReadData(Z, path);
00175 if(!fCrossSection[Z]) { return cs; }
00176 }
00177
00178 G4int idx = fNShells[Z]*6 - 4;
00179 if (gammaEnergy <= (fParam[Z])[idx-1]) { return cs; }
00180
00181 G4double x1 = 1.0/gammaEnergy;
00182 G4double x2 = x1*x1;
00183 G4double x3 = x2*x1;
00184
00185
00186 if(gammaEnergy >= (fParam[Z])[0]) {
00187 G4double x4 = x2*x2;
00188 cs = x1*((fParam[Z])[idx] + x1*(fParam[Z])[idx+1]
00189 + x2*(fParam[Z])[idx+2] + x3*(fParam[Z])[idx+3]
00190 + x4*(fParam[Z])[idx+4]);
00191
00192 } else if(gammaEnergy >= (fParam[Z])[1]) {
00193 cs = x3*(fCrossSection[Z])->Value(gammaEnergy);
00194
00195
00196 } else {
00197 cs = x3*(fCrossSectionLE[Z])->Value(gammaEnergy);
00198 }
00199 if (verboseLevel > 1) {
00200 G4cout << "LivermorePhotoElectricModel: E(keV)= " << gammaEnergy/keV
00201 << " Z= " << Z << " cross(barn)= " << cs/barn << G4endl;
00202 }
00203 return cs;
00204 }
00205
00206
00207
00208 void
00209 G4LivermorePhotoElectricModel::SampleSecondaries(
00210 std::vector<G4DynamicParticle*>* fvect,
00211 const G4MaterialCutsCouple* couple,
00212 const G4DynamicParticle* aDynamicGamma,
00213 G4double,
00214 G4double)
00215 {
00216 G4double gammaEnergy = aDynamicGamma->GetKineticEnergy();
00217 if (verboseLevel > 3) {
00218 G4cout << "G4LivermorePhotoElectricModel::SampleSecondaries() Egamma(keV)= "
00219 << gammaEnergy/keV << G4endl;
00220 }
00221
00222
00223 fParticleChange->SetProposedKineticEnergy(0.);
00224 fParticleChange->ProposeTrackStatus(fStopAndKill);
00225
00226
00227 G4ThreeVector photonDirection = aDynamicGamma->GetMomentumDirection();
00228
00229
00230
00231 const G4Element* elm = SelectRandomAtom(couple->GetMaterial(),theGamma,
00232 gammaEnergy);
00233 G4int Z = G4lrint(elm->GetZ());
00234
00235
00236
00237
00238
00239 if(Z >= maxZ) { Z = maxZ-1; }
00240
00241
00242 if(!fCrossSection[Z]) {
00243 char* path = getenv("G4LEDATA");
00244 ReadData(Z, path);
00245 if(!fCrossSection[Z]) {
00246 fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy);
00247 return;
00248 }
00249 }
00250
00251
00252 size_t shellIdx = 0;
00253 size_t nn = fNShellsUsed[Z];
00254
00255 if(nn > 1) {
00256 if(gammaEnergy >= (fParam[Z])[0]) {
00257 G4double x1 = 1.0/gammaEnergy;
00258 G4double x2 = x1*x1;
00259 G4double x3 = x2*x1;
00260 G4double x4 = x3*x1;
00261 G4int idx = nn*6 - 4;
00262
00263
00264 G4double cs0 = G4UniformRand()*((fParam[Z])[idx] + x1*(fParam[Z])[idx+1]
00265 + x2*(fParam[Z])[idx+2]
00266 + x3*(fParam[Z])[idx+3]
00267 + x4*(fParam[Z])[idx+4]);
00268 for(shellIdx=0; shellIdx<nn; ++shellIdx) {
00269 idx = shellIdx*6 + 2;
00270 if(gammaEnergy > (fParam[Z])[idx-1]) {
00271 G4double cs = (fParam[Z])[idx] + x1*(fParam[Z])[idx+1]
00272 + x2*(fParam[Z])[idx+2] + x3*(fParam[Z])[idx+3]
00273 + x4*(fParam[Z])[idx+4];
00274 if(cs >= cs0) { break; }
00275 }
00276 }
00277 if(shellIdx >= nn) { shellIdx = nn-1; }
00278
00279 } else {
00280
00281
00282
00283 G4double cs = G4UniformRand();
00284
00285 if(gammaEnergy >= (fParam[Z])[1]) {
00286 cs *= (fCrossSection[Z])->Value(gammaEnergy);
00287 } else {
00288 cs *= (fCrossSectionLE[Z])->Value(gammaEnergy);
00289 }
00290
00291 for(size_t j=0; j<nn; ++j) {
00292 shellIdx = (size_t)fShellCrossSection.GetComponentID(Z, j);
00293 if(gammaEnergy > (fParam[Z])[6*shellIdx+1]) {
00294 cs -= fShellCrossSection.GetValueForComponent(Z, j, gammaEnergy);
00295 }
00296 if(cs <= 0.0 || j+1 == nn) { break; }
00297 }
00298 }
00299 }
00300
00301 G4double bindingEnergy = (fParam[Z])[shellIdx*6 + 1];
00302
00303
00304
00305
00306
00307 const G4AtomicShell* shell = 0;
00308
00309
00310 if(fDeexcitationActive && shellIdx + 1 < nn) {
00311 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellIdx);
00312 shell = fAtomDeexcitation->GetAtomicShell(Z, as);
00313 }
00314
00315
00316
00317 if(gammaEnergy < bindingEnergy) {
00318 fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy);
00319 return;
00320 }
00321
00322
00323 G4double eKineticEnergy = gammaEnergy - bindingEnergy;
00324 G4double edep = bindingEnergy;
00325
00326
00327 G4ThreeVector electronDirection =
00328 GetAngularDistribution()->SampleDirection(aDynamicGamma,
00329 eKineticEnergy,
00330 shellIdx,
00331 couple->GetMaterial());
00332
00333
00334 G4DynamicParticle* electron = new G4DynamicParticle (theElectron,
00335 electronDirection,
00336 eKineticEnergy);
00337 fvect->push_back(electron);
00338
00339
00340 if(shell) {
00341 G4int index = couple->GetIndex();
00342 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
00343 size_t nbefore = fvect->size();
00344
00345 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
00346 size_t nafter = fvect->size();
00347 if(nafter > nbefore) {
00348 for (size_t j=nbefore; j<nafter; ++j) {
00349 edep -= ((*fvect)[j])->GetKineticEnergy();
00350 }
00351 }
00352 }
00353 }
00354
00355 if(edep > 0.0) {
00356 fParticleChange->ProposeLocalEnergyDeposit(edep);
00357 }
00358 }
00359
00360
00361
00362 void
00363 G4LivermorePhotoElectricModel::ReadData(G4int Z, const char* path)
00364 {
00365 if (verboseLevel > 1)
00366 {
00367 G4cout << "Calling ReadData() of G4LivermoreGammaConversionModel"
00368 << G4endl;
00369 }
00370
00371 if(fCrossSection[Z]) { return; }
00372
00373 const char* datadir = path;
00374
00375 if(!datadir)
00376 {
00377 datadir = getenv("G4LEDATA");
00378 if(!datadir)
00379 {
00380 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
00381 "em0006",FatalException,
00382 "Environment variable G4LEDATA not defined");
00383 return;
00384 }
00385 }
00386
00387
00388 fCrossSection[Z] = new G4LPhysicsFreeVector();
00389 fCrossSection[Z]->SetSpline(true);
00390
00391 std::ostringstream ost;
00392 ost << datadir << "/livermore/phot/pe-cs-" << Z <<".dat";
00393 std::ifstream fin(ost.str().c_str());
00394 if( !fin.is_open()) {
00395 G4ExceptionDescription ed;
00396 ed << "G4LivermorePhotoElectricModel data file <" << ost.str().c_str()
00397 << "> is not opened!" << G4endl;
00398 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
00399 "em0003",FatalException,
00400 ed,"G4LEDATA version should be G4EMLOW6.32 or later.");
00401 return;
00402 } else {
00403 if(verboseLevel > 3) { G4cout << "File " << ost.str().c_str()
00404 << " is opened by G4LivermorePhotoElectricModel" << G4endl;}
00405 fCrossSection[Z]->Retrieve(fin, true);
00406 fCrossSection[Z]->ScaleVector(MeV, barn);
00407 fin.close();
00408 }
00409
00410
00411 G4int n1 = 0;
00412 G4int n2 = 0;
00413 G4double x;
00414 std::ostringstream ost1;
00415 ost1 << datadir << "/livermore/phot/pe-" << Z <<".dat";
00416 std::ifstream fin1(ost1.str().c_str());
00417 if( !fin1.is_open()) {
00418 G4ExceptionDescription ed;
00419 ed << "G4LivermorePhotoElectricModel data file <" << ost1.str().c_str()
00420 << "> is not opened!" << G4endl;
00421 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
00422 "em0003",FatalException,
00423 ed,"G4LEDATA version should be G4EMLOW6.32 or later.");
00424 return;
00425 } else {
00426 if(verboseLevel > 3) {
00427 G4cout << "File " << ost1.str().c_str()
00428 << " is opened by G4LivermorePhotoElectricModel" << G4endl;
00429 }
00430 fin1 >> n1 >> n2 >> x;
00431 fNShells[Z] = n1;
00432 (fParam[Z]).reserve(6*n1+1);
00433 (fParam[Z]).push_back(x*MeV);
00434 for(G4int i=0; i<n1; ++i) {
00435 for(G4int j=0; j<6; ++j) {
00436 fin1 >> x;
00437 if(0 == j) { x *= MeV; }
00438 else { x *= barn; }
00439 (fParam[Z]).push_back(x);
00440 }
00441 }
00442 fin1.close();
00443 }
00444
00445 if(nShellLimit < n2) { n2 = nShellLimit; }
00446 fShellCrossSection.InitialiseForComponent(Z, n2);
00447 fNShellsUsed[Z] = n2;
00448
00449 if(1 < n2) {
00450 std::ostringstream ost2;
00451 ost2 << datadir << "/livermore/phot/pe-ss-cs-" << Z <<".dat";
00452 std::ifstream fin2(ost2.str().c_str());
00453 if( !fin2.is_open()) {
00454 G4ExceptionDescription ed;
00455 ed << "G4LivermorePhotoElectricModel data file <" << ost2.str().c_str()
00456 << "> is not opened!" << G4endl;
00457 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
00458 "em0003",FatalException,
00459 ed,"G4LEDATA version should be G4EMLOW6.32 or later.");
00460 return;
00461 } else {
00462 if(verboseLevel > 3) {
00463 G4cout << "File " << ost2.str().c_str()
00464 << " is opened by G4LivermorePhotoElectricModel" << G4endl;
00465 }
00466
00467 G4int n3, n4;
00468 G4double y;
00469 for(G4int i=0; i<n2; ++i) {
00470 fin2 >> x >> y >> n3 >> n4;
00471 G4LPhysicsFreeVector* v = new G4LPhysicsFreeVector(n3, x, y);
00472 for(G4int j=0; j<n3; ++j) {
00473 fin2 >> x >> y;
00474 v->PutValues(j, x*MeV, y*barn);
00475 }
00476 fShellCrossSection.AddComponent(Z, n4, v);
00477 }
00478 fin2.close();
00479 }
00480 }
00481
00482
00483 if(1 < fNShells[Z]) {
00484 fCrossSectionLE[Z] = new G4LPhysicsFreeVector();
00485
00486 std::ostringstream ost3;
00487 ost3 << datadir << "/livermore/phot/pe-le-cs-" << Z <<".dat";
00488 std::ifstream fin3(ost3.str().c_str());
00489 if( !fin3.is_open()) {
00490 G4ExceptionDescription ed;
00491 ed << "G4LivermorePhotoElectricModel data file <" << ost3.str().c_str()
00492 << "> is not opened!" << G4endl;
00493 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
00494 "em0003",FatalException,
00495 ed,"G4LEDATA version should be G4EMLOW6.32 or later.");
00496 return;
00497 } else {
00498 if(verboseLevel > 3) {
00499 G4cout << "File " << ost3.str().c_str()
00500 << " is opened by G4LivermorePhotoElectricModel" << G4endl;
00501 }
00502 fCrossSectionLE[Z]->Retrieve(fin3, true);
00503 fCrossSectionLE[Z]->ScaleVector(MeV, barn);
00504 fin3.close();
00505 }
00506 }
00507 }
00508
00509