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 #include "G4LivermoreNuclearGammaConversionModel.hh"
00032 #include "G4PhysicalConstants.hh"
00033 #include "G4SystemOfUnits.hh"
00034
00035
00036
00037 using namespace std;
00038
00039
00040
00041 G4LivermoreNuclearGammaConversionModel::G4LivermoreNuclearGammaConversionModel(const G4ParticleDefinition*,
00042 const G4String& nam)
00043 :G4VEmModel(nam),fParticleChange(0),smallEnergy(2.*MeV),
00044 isInitialised(false),
00045 crossSectionHandler(0),meanFreePathTable(0)
00046 {
00047 lowEnergyLimit = 2.0*electron_mass_c2;
00048 highEnergyLimit = 100 * GeV;
00049 SetHighEnergyLimit(highEnergyLimit);
00050
00051 verboseLevel= 0;
00052
00053
00054
00055
00056
00057
00058
00059 if(verboseLevel > 0) {
00060 G4cout << "Livermore Nuclear Gamma conversion is constructed " << G4endl
00061 << "Energy range: "
00062 << lowEnergyLimit / MeV << " MeV - "
00063 << highEnergyLimit / GeV << " GeV"
00064 << G4endl;
00065 }
00066 }
00067
00068
00069
00070 G4LivermoreNuclearGammaConversionModel::~G4LivermoreNuclearGammaConversionModel()
00071 {
00072 if (crossSectionHandler) delete crossSectionHandler;
00073 }
00074
00075
00076
00077 void
00078 G4LivermoreNuclearGammaConversionModel::Initialise(const G4ParticleDefinition*,
00079 const G4DataVector&)
00080 {
00081 if (verboseLevel > 3)
00082 G4cout << "Calling G4LivermoreNuclearGammaConversionModel::Initialise()" << G4endl;
00083
00084 if (crossSectionHandler)
00085 {
00086 crossSectionHandler->Clear();
00087 delete crossSectionHandler;
00088 }
00089
00090
00091
00092 crossSectionHandler = new G4CrossSectionHandler();
00093 crossSectionHandler->Initialise(0,lowEnergyLimit,100.*GeV,400);
00094 G4String crossSectionFile = "pairdata/pp-pair-cs-";
00095 crossSectionHandler->LoadData(crossSectionFile);
00096
00097
00098
00099 if (verboseLevel > 0) {
00100 G4cout << "Loaded cross section files for Livermore GammaConversion" << G4endl;
00101 G4cout << "To obtain the total cross section this should be used only " << G4endl
00102 << "in connection with G4ElectronGammaConversion " << G4endl;
00103 }
00104
00105 if (verboseLevel > 0) {
00106 G4cout << "Livermore Nuclear Gamma Conversion model is initialized " << G4endl
00107 << "Energy range: "
00108 << LowEnergyLimit() / MeV << " MeV - "
00109 << HighEnergyLimit() / GeV << " GeV"
00110 << G4endl;
00111 }
00112
00113 if(isInitialised) return;
00114 fParticleChange = GetParticleChangeForGamma();
00115 isInitialised = true;
00116 }
00117
00118
00119
00120 G4double
00121 G4LivermoreNuclearGammaConversionModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
00122 G4double GammaEnergy,
00123 G4double Z, G4double,
00124 G4double, G4double)
00125 {
00126 if (verboseLevel > 3) {
00127 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreNuclearGammaConversionModel"
00128 << G4endl;
00129 }
00130 if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0;
00131
00132 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
00133 return cs;
00134 }
00135
00136
00137
00138 void G4LivermoreNuclearGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00139 const G4MaterialCutsCouple* couple,
00140 const G4DynamicParticle* aDynamicGamma,
00141 G4double,
00142 G4double)
00143 {
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155 if (verboseLevel > 3)
00156 G4cout << "Calling SampleSecondaries() of G4LivermoreNuclearGammaConversionModel" << G4endl;
00157
00158 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
00159 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
00160
00161 G4double epsilon ;
00162 G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
00163
00164
00165 if (photonEnergy < smallEnergy )
00166 {
00167 epsilon = epsilon0Local + (0.5 - epsilon0Local) * G4UniformRand();
00168 }
00169 else
00170 {
00171
00172
00173 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
00174 const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
00175
00176 if (element == 0)
00177 {
00178 G4cout << "G4LivermoreNuclearGammaConversionModel::SampleSecondaries - element = 0"
00179 << G4endl;
00180 return;
00181 }
00182 G4IonisParamElm* ionisation = element->GetIonisation();
00183 if (ionisation == 0)
00184 {
00185 G4cout << "G4LivermoreNuclearGammaConversionModel::SampleSecondaries - ionisation = 0"
00186 << G4endl;
00187 return;
00188 }
00189
00190
00191 G4double fZ = 8. * (ionisation->GetlogZ3());
00192 if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb());
00193
00194
00195 G4double screenFactor = 136. * epsilon0Local / (element->GetIonisation()->GetZ3()) ;
00196 G4double screenMax = std::exp ((42.24 - fZ)/8.368) - 0.952 ;
00197 G4double screenMin = std::min(4.*screenFactor,screenMax) ;
00198
00199
00200 G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
00201 G4double epsilonMin = std::max(epsilon0Local,epsilon1);
00202 G4double epsilonRange = 0.5 - epsilonMin ;
00203
00204
00205 G4double screen;
00206 G4double gReject ;
00207
00208 G4double f10 = ScreenFunction1(screenMin) - fZ;
00209 G4double f20 = ScreenFunction2(screenMin) - fZ;
00210 G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
00211 G4double normF2 = std::max(1.5 * f20,0.);
00212
00213 do {
00214 if (normF1 / (normF1 + normF2) > G4UniformRand() )
00215 {
00216 epsilon = 0.5 - epsilonRange * std::pow(G4UniformRand(), 0.3333) ;
00217 screen = screenFactor / (epsilon * (1. - epsilon));
00218 gReject = (ScreenFunction1(screen) - fZ) / f10 ;
00219 }
00220 else
00221 {
00222 epsilon = epsilonMin + epsilonRange * G4UniformRand();
00223 screen = screenFactor / (epsilon * (1 - epsilon));
00224 gReject = (ScreenFunction2(screen) - fZ) / f20 ;
00225 }
00226 } while ( gReject < G4UniformRand() );
00227
00228 }
00229
00230
00231
00232 G4double electronTotEnergy;
00233 G4double positronTotEnergy;
00234
00235 if (G4int(2*G4UniformRand()))
00236 {
00237 electronTotEnergy = (1. - epsilon) * photonEnergy;
00238 positronTotEnergy = epsilon * photonEnergy;
00239 }
00240 else
00241 {
00242 positronTotEnergy = (1. - epsilon) * photonEnergy;
00243 electronTotEnergy = epsilon * photonEnergy;
00244 }
00245
00246
00247
00248
00249
00250 G4double u;
00251 const G4double a1 = 0.625;
00252 G4double a2 = 3. * a1;
00253
00254
00255
00256 if (0.25 > G4UniformRand())
00257 {
00258 u = - std::log(G4UniformRand() * G4UniformRand()) / a1 ;
00259 }
00260 else
00261 {
00262 u = - std::log(G4UniformRand() * G4UniformRand()) / a2 ;
00263 }
00264
00265 G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
00266 G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
00267 G4double phi = twopi * G4UniformRand();
00268
00269 G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
00270 G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
00271
00272
00273
00274
00275
00276
00277 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
00278
00279
00280
00281 G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
00282 electronDirection.rotateUz(photonDirection);
00283
00284 G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
00285 electronDirection,
00286 electronKineEnergy);
00287
00288
00289 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
00290
00291
00292
00293 G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
00294 positronDirection.rotateUz(photonDirection);
00295
00296
00297 G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
00298 positronDirection, positronKineEnergy);
00299
00300
00301 fvect->push_back(particle1);
00302 fvect->push_back(particle2);
00303
00304
00305 fParticleChange->SetProposedKineticEnergy(0.);
00306 fParticleChange->ProposeTrackStatus(fStopAndKill);
00307
00308 }
00309
00310
00311
00312 G4double G4LivermoreNuclearGammaConversionModel::ScreenFunction1(G4double screenVariable)
00313 {
00314
00315
00316 G4double value;
00317
00318 if (screenVariable > 1.)
00319 value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
00320 else
00321 value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
00322
00323 return value;
00324 }
00325
00326
00327
00328 G4double G4LivermoreNuclearGammaConversionModel::ScreenFunction2(G4double screenVariable)
00329 {
00330
00331
00332 G4double value;
00333
00334 if (screenVariable > 1.)
00335 value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
00336 else
00337 value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
00338
00339 return value;
00340 }
00341