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
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066 #include "G4LowEPComptonModel.hh"
00067 #include "G4PhysicalConstants.hh"
00068 #include "G4SystemOfUnits.hh"
00069 #include "G4Electron.hh"
00070 #include "G4ParticleChangeForGamma.hh"
00071 #include "G4LossTableManager.hh"
00072 #include "G4VAtomDeexcitation.hh"
00073 #include "G4AtomicShell.hh"
00074 #include "G4CrossSectionHandler.hh"
00075 #include "G4CompositeEMDataSet.hh"
00076 #include "G4LogLogInterpolation.hh"
00077 #include "G4Gamma.hh"
00078 #include "G4HadTmpUtil.hh"
00079
00080
00081 using namespace std;
00082
00083
00084
00085 G4LowEPComptonModel::G4LowEPComptonModel(const G4ParticleDefinition*,
00086 const G4String& nam)
00087 :G4VEmModel(nam),fParticleChange(0),isInitialised(false),
00088 scatterFunctionData(0),crossSectionHandler(0),fAtomDeexcitation(0)
00089 {
00090 lowEnergyLimit = 250 * eV;
00091 highEnergyLimit = 100 * GeV;
00092
00093 verboseLevel=0 ;
00094
00095
00096
00097
00098
00099
00100
00101 if( verboseLevel>0 ) {
00102 G4cout << "Low energy photon Compton model is constructed " << G4endl
00103 << "Energy range: "
00104 << lowEnergyLimit / eV << " eV - "
00105 << highEnergyLimit / GeV << " GeV"
00106 << G4endl;
00107 }
00108
00109
00110 SetDeexcitationFlag(true);
00111
00112 }
00113
00114
00115
00116 G4LowEPComptonModel::~G4LowEPComptonModel()
00117 {
00118 delete crossSectionHandler;
00119 delete scatterFunctionData;
00120 }
00121
00122
00123
00124 void G4LowEPComptonModel::Initialise(const G4ParticleDefinition* particle,
00125 const G4DataVector& cuts)
00126 {
00127 if (verboseLevel > 2) {
00128 G4cout << "Calling G4LowEPComptonModel::Initialise()" << G4endl;
00129 }
00130
00131 if (crossSectionHandler)
00132 {
00133 crossSectionHandler->Clear();
00134 delete crossSectionHandler;
00135 }
00136 delete scatterFunctionData;
00137
00138
00139 crossSectionHandler = new G4CrossSectionHandler;
00140 G4String crossSectionFile = "comp/ce-cs-";
00141 crossSectionHandler->LoadData(crossSectionFile);
00142
00143 G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
00144 G4String scatterFile = "comp/ce-sf-";
00145 scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
00146 scatterFunctionData->LoadData(scatterFile);
00147
00148
00149 shellData.SetOccupancyData();
00150 G4String file = "/doppler/shell-doppler";
00151 shellData.LoadData(file);
00152
00153 InitialiseElementSelectors(particle,cuts);
00154
00155 if (verboseLevel > 2) {
00156 G4cout << "Loaded cross section files for low energy photon Compton model" << G4endl;
00157 }
00158
00159 if(isInitialised) { return; }
00160 isInitialised = true;
00161
00162 fParticleChange = GetParticleChangeForGamma();
00163
00164 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
00165
00166 if( verboseLevel>0 ) {
00167 G4cout << "Low energy photon Compton model is initialized " << G4endl
00168 << "Energy range: "
00169 << LowEnergyLimit() / eV << " eV - "
00170 << HighEnergyLimit() / GeV << " GeV"
00171 << G4endl;
00172 }
00173 }
00174
00175
00176
00177 G4double G4LowEPComptonModel::ComputeCrossSectionPerAtom(
00178 const G4ParticleDefinition*,
00179 G4double GammaEnergy,
00180 G4double Z, G4double,
00181 G4double, G4double)
00182 {
00183 if (verboseLevel > 3) {
00184 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LowEPComptonModel" << G4endl;
00185 }
00186 if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) { return 0.0; }
00187
00188 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
00189 return cs;
00190 }
00191
00192
00193
00194
00195
00196
00197
00198
00199 void G4LowEPComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00200 const G4MaterialCutsCouple* couple,
00201 const G4DynamicParticle* aDynamicGamma,
00202 G4double, G4double)
00203 {
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy()/MeV;
00220
00221 if (verboseLevel > 3) {
00222 G4cout << "G4LowEPComptonModel::SampleSecondaries() E(MeV)= "
00223 << photonEnergy0/MeV << " in " << couple->GetMaterial()->GetName()
00224 << G4endl;
00225 }
00226
00227
00228 if (photonEnergy0 <= lowEnergyLimit)
00229 {
00230 fParticleChange->ProposeTrackStatus(fStopAndKill);
00231 fParticleChange->SetProposedKineticEnergy(0.);
00232 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
00233 return ;
00234 }
00235
00236 G4double e0m = photonEnergy0 / electron_mass_c2 ;
00237 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
00238
00239
00240 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
00241 const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
00242 G4int Z = (G4int)elm->GetZ();
00243
00244 G4double LowEPCepsilon0 = 1. / (1. + 2. * e0m);
00245 G4double LowEPCepsilon0Sq = LowEPCepsilon0 * LowEPCepsilon0;
00246 G4double alpha1 = -std::log(LowEPCepsilon0);
00247 G4double alpha2 = 0.5 * (1. - LowEPCepsilon0Sq);
00248
00249 G4double wlPhoton = h_Planck*c_light/photonEnergy0;
00250
00251
00252 G4double LowEPCepsilon;
00253 G4double LowEPCepsilonSq;
00254 G4double oneCosT;
00255 G4double sinT2;
00256 G4double gReject;
00257
00258 do
00259 {
00260 if ( alpha1/(alpha1+alpha2) > G4UniformRand())
00261 {
00262 LowEPCepsilon = std::exp(-alpha1 * G4UniformRand());
00263 LowEPCepsilonSq = LowEPCepsilon * LowEPCepsilon;
00264 }
00265 else
00266 {
00267 LowEPCepsilonSq = LowEPCepsilon0Sq + (1. - LowEPCepsilon0Sq) * G4UniformRand();
00268 LowEPCepsilon = std::sqrt(LowEPCepsilonSq);
00269 }
00270
00271 oneCosT = (1. - LowEPCepsilon) / ( LowEPCepsilon * e0m);
00272 sinT2 = oneCosT * (2. - oneCosT);
00273 G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm);
00274 G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
00275 gReject = (1. - LowEPCepsilon * sinT2 / (1. + LowEPCepsilonSq)) * scatteringFunction;
00276
00277 } while(gReject < G4UniformRand()*Z);
00278
00279 G4double cosTheta = 1. - oneCosT;
00280 G4double sinTheta = std::sqrt(sinT2);
00281 G4double phi = twopi * G4UniformRand();
00282 G4double dirx = sinTheta * std::cos(phi);
00283 G4double diry = sinTheta * std::sin(phi);
00284 G4double dirz = cosTheta ;
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294 G4double vel_c = 299792458;
00295 G4double momentum_au_to_nat = (pi/2.0)*1.992851740*std::pow(10.,-24.);
00296 G4double e_mass_kg = 9.10938188 * std::pow(10.,-31.);
00297
00298 G4int maxDopplerIterations = 1000;
00299 G4double bindingE = 0.;
00300 G4double pEIncident = photonEnergy0 ;
00301 G4double pERecoil = -1.;
00302 G4double eERecoil = -1.;
00303 G4double e_alpha =0.;
00304 G4double e_beta = 0.;
00305
00306 G4double CE_emission_flag = 0.;
00307 G4double ePAU = -1;
00308
00309 G4int shellIdx = 0;
00310 G4double u_temp = 0;
00311 G4double cosPhiE =0;
00312 G4double sinThetaE =0;
00313 G4double cosThetaE =0;
00314 G4int iteration = 0;
00315 do{
00316
00317
00318
00319
00320
00321
00322 do
00323 {
00324 iteration++;
00325
00326
00327
00328
00329
00330
00331
00332
00333 shellIdx = shellData.SelectRandomShell(Z);
00334 bindingE = shellData.BindingEnergy(Z,shellIdx)/MeV;
00335
00336
00337
00338
00339 ePAU = profileData.RandomSelectMomentum(Z,shellIdx);
00340
00341
00342
00343 G4double ePSI = ePAU * momentum_au_to_nat;
00344
00345
00346
00347 u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / ((e_mass_kg*e_mass_kg)*(vel_c*vel_c)+(ePSI*ePSI)))/vel_c;
00348
00349
00350
00351 e_alpha = pi*G4UniformRand();
00352 e_beta = twopi*G4UniformRand();
00353
00354
00355
00356 G4double eEIncident = electron_mass_c2 / sqrt( 1 - (u_temp*u_temp));
00357 G4double systemE = eEIncident + pEIncident;
00358
00359
00360 G4double gamma_temp = 1.0 / sqrt( 1 - (u_temp*u_temp));
00361 G4double numerator = gamma_temp*electron_mass_c2*(1 - u_temp * std::cos(e_alpha));
00362 G4double subdenom1 = u_temp*cosTheta*std::cos(e_alpha);
00363 G4double subdenom2 = u_temp*sinTheta*std::sin(e_alpha)*std::cos(e_beta);
00364 G4double denominator = (1.0 - cosTheta) + (gamma_temp*electron_mass_c2*(1 - subdenom1 - subdenom2) / pEIncident);
00365 pERecoil = (numerator/denominator);
00366 eERecoil = systemE - pERecoil;
00367 CE_emission_flag = pEIncident - pERecoil;
00368 } while ( (iteration <= maxDopplerIterations) && (CE_emission_flag < bindingE));
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382 G4double a_temp = eERecoil / electron_mass_c2;
00383 G4double u_p_temp = sqrt(1 - (1 / (a_temp*a_temp)));
00384
00385
00386
00387 G4double sinAlpha = std::sin(e_alpha);
00388 G4double cosAlpha = std::cos(e_alpha);
00389 G4double sinBeta = std::sin(e_beta);
00390 G4double cosBeta = std::cos(e_beta);
00391
00392 G4double gamma = 1.0 / sqrt(1 - (u_temp*u_temp));
00393 G4double gamma_p = 1.0 / sqrt(1 - (u_p_temp*u_p_temp));
00394
00395 G4double var_A = pERecoil*u_p_temp*sinTheta;
00396 G4double var_B = u_p_temp* (pERecoil*cosTheta-pEIncident);
00397 G4double var_C = (pERecoil-pEIncident) - ( (pERecoil*pEIncident) / (gamma_p*electron_mass_c2))*(1 - cosTheta);
00398
00399 G4double var_D1 = gamma*electron_mass_c2*pERecoil;
00400 G4double var_D2 = (1 - (u_temp*cosTheta*cosAlpha) - (u_temp*sinTheta*cosBeta*sinAlpha));
00401 G4double var_D3 = ((electron_mass_c2*electron_mass_c2)*(gamma*gamma_p - 1)) - (gamma_p*electron_mass_c2*pERecoil);
00402 G4double var_D = var_D1*var_D2 + var_D3;
00403
00404 G4double var_E1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosAlpha);
00405 G4double var_E2 = gamma_p*electron_mass_c2*pERecoil*u_p_temp*cosTheta;
00406 G4double var_E = var_E1 - var_E2;
00407
00408 G4double var_F1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosBeta*sinAlpha);
00409 G4double var_F2 = (gamma_p*electron_mass_c2*pERecoil*u_p_temp*sinTheta);
00410 G4double var_F = var_F1 - var_F2;
00411
00412 G4double var_G = (gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*sinBeta*sinAlpha;
00413
00414
00415
00416
00417 G4double var_W1 = (var_F*var_B - var_E*var_A)*(var_F*var_B - var_E*var_A);
00418 G4double var_W2 = (var_G*var_G)*(var_A*var_A) + (var_G*var_G)*(var_B*var_B);
00419 G4double var_W = var_W1 + var_W2;
00420
00421 G4double var_Y = 2.0*(((var_A*var_D-var_F*var_C)*(var_F*var_B-var_E*var_A)) - ((var_G*var_G)*var_B*var_C));
00422
00423 G4double var_Z1 = (var_A*var_D - var_F*var_C)*(var_A*var_D - var_F*var_C);
00424 G4double var_Z2 = (var_G*var_G)*(var_C*var_C) - (var_G*var_G)*(var_A*var_A);
00425 G4double var_Z = var_Z1 + var_Z2;
00426
00427 G4double diff = ((var_Y*var_Y)-4*var_W*var_Z);
00428
00429
00430
00431 if (diff < 0.0)
00432 {
00433
00434 G4double funorder = 0.0;
00435 G4double sf = 6.0;
00436
00437 G4double diff1 = var_Y*var_Y;
00438 funorder = abs(G4lrint(std::log10(static_cast<double>(abs(diff1)))+std::numeric_limits<double>::epsilon()))+sf;
00439 diff1 = G4lrint(diff1*std::pow(10.0,funorder))/(1.0*std::pow(10.0,funorder));
00440
00441 G4double diff2 = 4*var_W*var_Z;
00442 funorder = abs(G4lrint(std::log10(static_cast<double>(abs(diff2)))+std::numeric_limits<double>::epsilon()))+sf;
00443 diff2 = G4lrint(diff2*std::pow(10.0,funorder))/(1.0*std::pow(10.0,funorder));
00444
00445 diff = diff1 -diff2;
00446
00447 }
00448
00449
00450
00451
00452 G4double X_p = (-var_Y + sqrt (diff))/(2*var_W);
00453 G4double X_m = (-var_Y - sqrt (diff))/(2*var_W);
00454
00455
00456
00457 G4double ThetaE = 0.;
00458 G4double sol_select = G4UniformRand();
00459
00460 if (sol_select < 0.5)
00461 {
00462 ThetaE = std::acos(X_p);
00463 }
00464 if (sol_select > 0.5)
00465 {
00466 ThetaE = std::acos(X_m);
00467 }
00468
00469 cosThetaE = std::cos(ThetaE);
00470 sinThetaE = std::sin(ThetaE);
00471 G4double Theta = std::acos(cosTheta);
00472
00473
00474 G4double iSinThetaE = std::sqrt(1+std::tan((pi/2.0)-ThetaE)*std::tan((pi/2.0)-ThetaE));
00475 G4double iSinTheta = std::sqrt(1+std::tan((pi/2.0)-Theta)*std::tan((pi/2.0)-Theta));
00476 G4double ivar_A = iSinTheta/ (pERecoil*u_p_temp);
00477
00478 cosPhiE = (var_C - var_B*cosThetaE)*(ivar_A*iSinThetaE);
00479
00480
00481
00482 if (abs(cosPhiE) > 1.0)
00483 {
00484 G4double funorder = 0.0;
00485 G4double sf = 4.0;
00486 funorder = abs(G4lrint(std::log10(static_cast<double>(abs(cosPhiE)))+std::numeric_limits<double>::epsilon()))+sf;
00487 cosPhiE = G4lrint(cosPhiE*std::pow(10.0,funorder))/(1.0*std::pow(10.0,funorder));
00488
00489 }
00490
00491
00492
00493
00494 } while ( (iteration <= maxDopplerIterations) && (abs(cosPhiE) > 1));
00495
00496
00497
00498
00499 if (iteration >= maxDopplerIterations)
00500 {
00501 pERecoil = photonEnergy0 ;
00502 bindingE = 0.;
00503 dirx=0.0;
00504 diry=0.0;
00505 dirz=1.0;
00506 }
00507
00508
00509
00510 G4ThreeVector photonDirection1(dirx,diry,dirz);
00511 photonDirection1.rotateUz(photonDirection0);
00512 fParticleChange->ProposeMomentumDirection(photonDirection1) ;
00513
00514 if (pERecoil > 0.)
00515 {
00516 fParticleChange->SetProposedKineticEnergy(pERecoil) ;
00517
00518
00519 G4double PhiE = std::acos(cosPhiE);
00520 G4double eDirX = sinThetaE * std::cos(phi+PhiE);
00521 G4double eDirY = sinThetaE * std::sin(phi+PhiE);
00522 G4double eDirZ = cosThetaE;
00523
00524 G4double eKineticEnergy = pEIncident - pERecoil - bindingE;
00525
00526 G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
00527 eDirection.rotateUz(photonDirection0);
00528 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),
00529 eDirection,eKineticEnergy) ;
00530 fvect->push_back(dp);
00531
00532 }
00533 else
00534 {
00535 fParticleChange->SetProposedKineticEnergy(0.);
00536 fParticleChange->ProposeTrackStatus(fStopAndKill);
00537 }
00538
00539
00540
00541
00542
00543
00544 if(fAtomDeexcitation && iteration < maxDopplerIterations) {
00545 G4int index = couple->GetIndex();
00546 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
00547 size_t nbefore = fvect->size();
00548 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellIdx);
00549 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
00550 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
00551 size_t nafter = fvect->size();
00552 if(nafter > nbefore) {
00553 for (size_t i=nbefore; i<nafter; ++i) {
00554 bindingE -= ((*fvect)[i])->GetKineticEnergy();
00555 }
00556 }
00557 }
00558 }
00559 if(bindingE < 0.0) { bindingE = 0.0; }
00560 fParticleChange->ProposeLocalEnergyDeposit(bindingE);
00561
00562 }
00563