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 "G4LivermoreComptonModifiedModel.hh"
00045 #include "G4PhysicalConstants.hh"
00046 #include "G4SystemOfUnits.hh"
00047 #include "G4Electron.hh"
00048 #include "G4ParticleChangeForGamma.hh"
00049 #include "G4LossTableManager.hh"
00050 #include "G4VAtomDeexcitation.hh"
00051 #include "G4AtomicShell.hh"
00052 #include "G4CrossSectionHandler.hh"
00053 #include "G4CompositeEMDataSet.hh"
00054 #include "G4LogLogInterpolation.hh"
00055 #include "G4Gamma.hh"
00056
00057
00058
00059 using namespace std;
00060
00061
00062
00063 G4LivermoreComptonModifiedModel::G4LivermoreComptonModifiedModel(const G4ParticleDefinition*,
00064 const G4String& nam)
00065 :G4VEmModel(nam),fParticleChange(0),isInitialised(false),
00066 scatterFunctionData(0),
00067 crossSectionHandler(0),fAtomDeexcitation(0)
00068 {
00069 lowEnergyLimit = 250 * eV;
00070 highEnergyLimit = 100 * GeV;
00071
00072 verboseLevel=0 ;
00073
00074
00075
00076
00077
00078
00079
00080 if( verboseLevel>0 ) {
00081 G4cout << "Livermore Modified Compton model is constructed " << G4endl
00082 << "Energy range: "
00083 << lowEnergyLimit / eV << " eV - "
00084 << highEnergyLimit / GeV << " GeV"
00085 << G4endl;
00086 }
00087
00088
00089 SetDeexcitationFlag(true);
00090
00091 }
00092
00093
00094
00095 G4LivermoreComptonModifiedModel::~G4LivermoreComptonModifiedModel()
00096 {
00097 delete crossSectionHandler;
00098 delete scatterFunctionData;
00099 }
00100
00101
00102
00103 void G4LivermoreComptonModifiedModel::Initialise(const G4ParticleDefinition* particle,
00104 const G4DataVector& cuts)
00105 {
00106 if (verboseLevel > 2) {
00107 G4cout << "Calling G4LivermoreComptonModifiedModel::Initialise()" << G4endl;
00108 }
00109
00110 if (crossSectionHandler)
00111 {
00112 crossSectionHandler->Clear();
00113 delete crossSectionHandler;
00114 }
00115 delete scatterFunctionData;
00116
00117
00118 crossSectionHandler = new G4CrossSectionHandler;
00119 G4String crossSectionFile = "comp/ce-cs-";
00120 crossSectionHandler->LoadData(crossSectionFile);
00121
00122 G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
00123 G4String scatterFile = "comp/ce-sf-";
00124 scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
00125 scatterFunctionData->LoadData(scatterFile);
00126
00127
00128 shellData.SetOccupancyData();
00129 G4String file = "/doppler/shell-doppler";
00130 shellData.LoadData(file);
00131
00132 InitialiseElementSelectors(particle,cuts);
00133
00134 if (verboseLevel > 2) {
00135 G4cout << "Loaded cross section files for Livermore Modified Compton model" << G4endl;
00136 }
00137
00138 if(isInitialised) { return; }
00139 isInitialised = true;
00140
00141 fParticleChange = GetParticleChangeForGamma();
00142
00143 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
00144
00145 if( verboseLevel>0 ) {
00146 G4cout << "Livermore Compton model is initialized " << G4endl
00147 << "Energy range: "
00148 << LowEnergyLimit() / eV << " eV - "
00149 << HighEnergyLimit() / GeV << " GeV"
00150 << G4endl;
00151 }
00152 }
00153
00154
00155
00156 G4double G4LivermoreComptonModifiedModel::ComputeCrossSectionPerAtom(
00157 const G4ParticleDefinition*,
00158 G4double GammaEnergy,
00159 G4double Z, G4double,
00160 G4double, G4double)
00161 {
00162 if (verboseLevel > 3) {
00163 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreComptonModifiedModel" << G4endl;
00164 }
00165 if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) { return 0.0; }
00166
00167 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
00168 return cs;
00169 }
00170
00171
00172
00173 void G4LivermoreComptonModifiedModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00174 const G4MaterialCutsCouple* couple,
00175 const G4DynamicParticle* aDynamicGamma,
00176 G4double, G4double)
00177 {
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
00194
00195 if (verboseLevel > 3) {
00196 G4cout << "G4LivermoreComptonModifiedModel::SampleSecondaries() E(MeV)= "
00197 << photonEnergy0/MeV << " in " << couple->GetMaterial()->GetName()
00198 << G4endl;
00199 }
00200
00201
00202 if (photonEnergy0 <= lowEnergyLimit)
00203 {
00204 fParticleChange->ProposeTrackStatus(fStopAndKill);
00205 fParticleChange->SetProposedKineticEnergy(0.);
00206 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
00207 return ;
00208 }
00209
00210 G4double e0m = photonEnergy0 / electron_mass_c2 ;
00211 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
00212
00213
00214 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
00215 const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
00216 G4int Z = (G4int)elm->GetZ();
00217
00218 G4double epsilon0Local = 1. / (1. + 2. * e0m);
00219 G4double epsilon0Sq = epsilon0Local * epsilon0Local;
00220 G4double alpha1 = -std::log(epsilon0Local);
00221 G4double alpha2 = 0.5 * (1. - epsilon0Sq);
00222
00223 G4double wlPhoton = h_Planck*c_light/photonEnergy0;
00224
00225
00226 G4double epsilon;
00227 G4double epsilonSq;
00228 G4double oneCosT;
00229 G4double sinT2;
00230 G4double gReject;
00231
00232 do
00233 {
00234 if ( alpha1/(alpha1+alpha2) > G4UniformRand())
00235 {
00236
00237 epsilon = std::exp(-alpha1 * G4UniformRand());
00238 epsilonSq = epsilon * epsilon;
00239 }
00240 else
00241 {
00242 epsilonSq = epsilon0Sq + (1. - epsilon0Sq) * G4UniformRand();
00243 epsilon = std::sqrt(epsilonSq);
00244 }
00245
00246 oneCosT = (1. - epsilon) / ( epsilon * e0m);
00247 sinT2 = oneCosT * (2. - oneCosT);
00248 G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm);
00249 G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
00250 gReject = (1. - epsilon * sinT2 / (1. + epsilonSq)) * scatteringFunction;
00251
00252 } while(gReject < G4UniformRand()*Z);
00253
00254 G4double cosTheta = 1. - oneCosT;
00255 G4double sinTheta = std::sqrt (sinT2);
00256 G4double phi = twopi * G4UniformRand() ;
00257 G4double dirx = sinTheta * std::cos(phi);
00258 G4double diry = sinTheta * std::sin(phi);
00259 G4double dirz = cosTheta ;
00260
00261
00262
00263
00264
00265
00266
00267 G4int maxDopplerIterations = 1000;
00268 G4double bindingE = 0.;
00269 G4double photonEoriginal = epsilon * photonEnergy0;
00270 G4double photonE = -1.;
00271 G4int iteration = 0;
00272 G4double systemE = 0;
00273 G4double ePAU = -1;
00274 G4int shellIdx = 0;
00275 G4double vel_c = 299792458;
00276 G4double momentum_au_to_nat = 1.992851740*std::pow(10.,-24.);
00277 G4double e_mass_kg = 9.10938188 * std::pow(10.,-31.);
00278 G4double eMax = -1;
00279 G4double Alpha=0;
00280 do
00281 {
00282 ++iteration;
00283
00284 shellIdx = shellData.SelectRandomShell(Z);
00285 bindingE = shellData.BindingEnergy(Z,shellIdx);
00286
00287
00288
00289
00290
00291 G4double pSample = profileData.RandomSelectMomentum(Z,shellIdx);
00292
00293
00294
00295
00296
00297
00298
00299
00300 do {
00301 Alpha = G4UniformRand()*pi/2.0;
00302 } while(Alpha >= (pi/2.0));
00303
00304 ePAU = pSample / std::cos(Alpha);
00305
00306
00307
00308 G4double ePSI = ePAU * momentum_au_to_nat;
00309 G4double u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / ((e_mass_kg*e_mass_kg)*(vel_c*vel_c)+(ePSI*ePSI)))/vel_c;
00310 G4double eEIncident = electron_mass_c2 / sqrt( 1 - (u_temp*u_temp));
00311
00312
00313 systemE = eEIncident+photonEnergy0;
00314
00315 eMax = systemE - bindingE - electron_mass_c2;
00316 G4double pDoppler = pSample * fine_structure_const;
00317 G4double pDoppler2 = pDoppler * pDoppler;
00318 G4double var2 = 1. + oneCosT * e0m;
00319 G4double var3 = var2*var2 - pDoppler2;
00320 G4double var4 = var2 - pDoppler2 * cosTheta;
00321 G4double var = var4*var4 - var3 + pDoppler2 * var3;
00322 if (var > 0.)
00323 {
00324 G4double varSqrt = std::sqrt(var);
00325 G4double scale = photonEnergy0 / var3;
00326
00327 if (G4UniformRand() < 0.5) { photonE = (var4 - varSqrt) * scale; }
00328 else { photonE = (var4 + varSqrt) * scale; }
00329 }
00330 else
00331 {
00332 photonE = -1.;
00333 }
00334 } while ( iteration <= maxDopplerIterations &&
00335 (photonE < 0. || photonE > eMax ) );
00336
00337
00338
00339 G4double eKineticEnergy = systemE - photonE - bindingE - electron_mass_c2;
00340
00341
00342 G4double eDirX = 0.0;
00343 G4double eDirY = 0.0;
00344 G4double eDirZ = 1.0;
00345
00346 if(eKineticEnergy < 0.0) {
00347 G4cout << "Error, kinetic energy of electron less than zero" << G4endl;
00348 }
00349
00350 else{
00351
00352
00353
00354 G4double E_num = photonEnergy0 - photonE*cosTheta;
00355 G4double E_dom = sqrt(photonEnergy0*photonEnergy0 + photonE*photonE -2*photonEnergy0*photonE*cosTheta);
00356 G4double cosThetaE = E_num / E_dom;
00357 G4double sinThetaE = -sqrt((1. - cosThetaE) * (1. + cosThetaE));
00358
00359 eDirX = sinThetaE * std::cos(phi);
00360 eDirY = sinThetaE * std::sin(phi);
00361 eDirZ = cosThetaE;
00362
00363 G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
00364 eDirection.rotateUz(photonDirection0);
00365 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),
00366 eDirection,eKineticEnergy) ;
00367 fvect->push_back(dp);
00368 }
00369
00370
00371
00372
00373 if (iteration >= maxDopplerIterations)
00374 {
00375 photonE = photonEoriginal;
00376 bindingE = 0.;
00377 }
00378
00379
00380
00381 G4ThreeVector photonDirection1(dirx,diry,dirz);
00382 photonDirection1.rotateUz(photonDirection0);
00383 fParticleChange->ProposeMomentumDirection(photonDirection1) ;
00384
00385 G4double photonEnergy1 = photonE;
00386
00387 if (photonEnergy1 > 0.)
00388 {
00389 fParticleChange->SetProposedKineticEnergy(photonEnergy1) ;
00390
00391 if (iteration < maxDopplerIterations)
00392 {
00393 G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
00394 eDirection.rotateUz(photonDirection0);
00395 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),
00396 eDirection,eKineticEnergy) ;
00397 fvect->push_back(dp);
00398 }
00399 }
00400 else
00401 {
00402 photonEnergy1 = 0.;
00403 fParticleChange->SetProposedKineticEnergy(0.) ;
00404 fParticleChange->ProposeTrackStatus(fStopAndKill);
00405 }
00406
00407
00408
00409 if(fAtomDeexcitation && iteration < maxDopplerIterations) {
00410 G4int index = couple->GetIndex();
00411 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
00412 size_t nbefore = fvect->size();
00413 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellIdx);
00414 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
00415 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
00416 size_t nafter = fvect->size();
00417 if(nafter > nbefore) {
00418 for (size_t i=nbefore; i<nafter; ++i) {
00419 bindingE -= ((*fvect)[i])->GetKineticEnergy();
00420 }
00421 }
00422 }
00423 }
00424 if(bindingE < 0.0) { bindingE = 0.0; }
00425 fParticleChange->ProposeLocalEnergyDeposit(bindingE);
00426 }
00427