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 "G4LivermoreComptonModel.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 G4LivermoreComptonModel::G4LivermoreComptonModel(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 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 G4LivermoreComptonModel::~G4LivermoreComptonModel()
00096 {
00097 delete crossSectionHandler;
00098 delete scatterFunctionData;
00099 }
00100
00101
00102
00103 void G4LivermoreComptonModel::Initialise(const G4ParticleDefinition* particle,
00104 const G4DataVector& cuts)
00105 {
00106 if (verboseLevel > 2) {
00107 G4cout << "Calling G4LivermoreComptonModel::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 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 G4LivermoreComptonModel::ComputeCrossSectionPerAtom(
00157 const G4ParticleDefinition*,
00158 G4double GammaEnergy,
00159 G4double Z, G4double,
00160 G4double, G4double)
00161 {
00162 if (verboseLevel > 3) {
00163 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreComptonModel" << 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 G4LivermoreComptonModel::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 << "G4LivermoreComptonModel::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 eMax = photonEnergy0;
00273
00274 G4int shellIdx = 0;
00275
00276 do
00277 {
00278 ++iteration;
00279
00280 shellIdx = shellData.SelectRandomShell(Z);
00281 bindingE = shellData.BindingEnergy(Z,shellIdx);
00282
00283 eMax = photonEnergy0 - bindingE;
00284
00285
00286
00287 G4double pSample = profileData.RandomSelectMomentum(Z,shellIdx);
00288
00289 G4double pDoppler = pSample * fine_structure_const;
00290 G4double pDoppler2 = pDoppler * pDoppler;
00291 G4double var2 = 1. + oneCosT * e0m;
00292 G4double var3 = var2*var2 - pDoppler2;
00293 G4double var4 = var2 - pDoppler2 * cosTheta;
00294 G4double var = var4*var4 - var3 + pDoppler2 * var3;
00295 if (var > 0.)
00296 {
00297 G4double varSqrt = std::sqrt(var);
00298 G4double scale = photonEnergy0 / var3;
00299
00300 if (G4UniformRand() < 0.5) { photonE = (var4 - varSqrt) * scale; }
00301 else { photonE = (var4 + varSqrt) * scale; }
00302 }
00303 else
00304 {
00305 photonE = -1.;
00306 }
00307 } while ( iteration <= maxDopplerIterations &&
00308
00309
00310
00311
00312 (photonE > eMax ) );
00313
00314
00315
00316
00317 if (iteration >= maxDopplerIterations)
00318 {
00319 photonE = photonEoriginal;
00320 bindingE = 0.;
00321 }
00322
00323
00324
00325 G4ThreeVector photonDirection1(dirx,diry,dirz);
00326 photonDirection1.rotateUz(photonDirection0);
00327 fParticleChange->ProposeMomentumDirection(photonDirection1) ;
00328
00329 G4double photonEnergy1 = photonE;
00330
00331 if (photonEnergy1 > 0.)
00332 {
00333 fParticleChange->SetProposedKineticEnergy(photonEnergy1) ;
00334 }
00335 else
00336 {
00337 photonEnergy1 = 0.;
00338 fParticleChange->SetProposedKineticEnergy(0.) ;
00339 fParticleChange->ProposeTrackStatus(fStopAndKill);
00340 }
00341
00342
00343 G4double eKineticEnergy = photonEnergy0 - photonEnergy1 - bindingE;
00344
00345
00346 if(eKineticEnergy < 0.0) {
00347 bindingE = photonEnergy0 - photonEnergy1;
00348
00349 } else {
00350 G4double eTotalEnergy = eKineticEnergy + electron_mass_c2;
00351
00352 G4double electronE = photonEnergy0 * (1. - epsilon) + electron_mass_c2;
00353 G4double electronP2 = electronE*electronE - electron_mass_c2*electron_mass_c2;
00354 G4double sinThetaE = -1.;
00355 G4double cosThetaE = 0.;
00356 if (electronP2 > 0.)
00357 {
00358 cosThetaE = (eTotalEnergy + photonEnergy1 )* (1. - epsilon) / std::sqrt(electronP2);
00359 sinThetaE = -sqrt((1. - cosThetaE) * (1. + cosThetaE));
00360 }
00361
00362 G4double eDirX = sinThetaE * std::cos(phi);
00363 G4double eDirY = sinThetaE * std::sin(phi);
00364 G4double eDirZ = cosThetaE;
00365
00366 G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
00367 eDirection.rotateUz(photonDirection0);
00368 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),
00369 eDirection,eKineticEnergy) ;
00370 fvect->push_back(dp);
00371 }
00372
00373
00374
00375 if(fAtomDeexcitation && iteration < maxDopplerIterations) {
00376 G4int index = couple->GetIndex();
00377 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
00378 size_t nbefore = fvect->size();
00379 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellIdx);
00380 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
00381 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
00382 size_t nafter = fvect->size();
00383 if(nafter > nbefore) {
00384 for (size_t i=nbefore; i<nafter; ++i) {
00385 bindingE -= ((*fvect)[i])->GetKineticEnergy();
00386 }
00387 }
00388 }
00389 }
00390 if(bindingE < 0.0) { bindingE = 0.0; }
00391 fParticleChange->ProposeLocalEnergyDeposit(bindingE);
00392 }
00393