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 #include "G4LivermoreGammaConversionModel.hh"
00031 #include "G4PhysicalConstants.hh"
00032 #include "G4SystemOfUnits.hh"
00033
00034
00035
00036 using namespace std;
00037
00038
00039
00040 G4LivermoreGammaConversionModel::G4LivermoreGammaConversionModel(const G4ParticleDefinition*,
00041 const G4String& nam)
00042 :G4VEmModel(nam),smallEnergy(2.*MeV),isInitialised(false),maxZ(99)
00043 {
00044 fParticleChange = 0;
00045
00046 lowEnergyLimit = 2.0*electron_mass_c2;
00047 data.resize(maxZ+1,0);
00048
00049 verboseLevel= 0;
00050
00051
00052
00053
00054
00055 if(verboseLevel > 0)
00056 {
00057 G4cout << "G4LivermoreGammaConversionModel is constructed " << G4endl;
00058 }
00059 }
00060
00061
00062
00063 G4LivermoreGammaConversionModel::~G4LivermoreGammaConversionModel()
00064 {
00065 for(G4int i=0; i<=maxZ; ++i) { delete data[i]; }
00066 }
00067
00068
00069
00070 void
00071 G4LivermoreGammaConversionModel::Initialise(const G4ParticleDefinition* particle,
00072 const G4DataVector& cuts)
00073 {
00074 if (verboseLevel > 1)
00075 {
00076 G4cout << "Calling Initialise() of G4LivermoreGammaConversionModel." << G4endl
00077 << "Energy range: "
00078 << LowEnergyLimit() / MeV << " MeV - "
00079 << HighEnergyLimit() / GeV << " GeV"
00080 << G4endl;
00081 }
00082
00083
00084
00085 InitialiseElementSelectors(particle, cuts);
00086
00087
00088
00089 char* path = getenv("G4LEDATA");
00090
00091 G4ProductionCutsTable* theCoupleTable =
00092 G4ProductionCutsTable::GetProductionCutsTable();
00093 G4int numOfCouples = theCoupleTable->GetTableSize();
00094
00095 for(G4int i=0; i<numOfCouples; ++i)
00096 {
00097 const G4Material* material =
00098 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
00099 const G4ElementVector* theElementVector = material->GetElementVector();
00100 G4int nelm = material->GetNumberOfElements();
00101
00102 for (G4int j=0; j<nelm; ++j)
00103 {
00104
00105 G4int Z = (G4int)(*theElementVector)[j]->GetZ();
00106 if(Z < 1) { Z = 1; }
00107 else if(Z > maxZ) { Z = maxZ; }
00108 if(!data[Z]) { ReadData(Z, path); }
00109 }
00110 }
00111
00112
00113 if(isInitialised) { return; }
00114 fParticleChange = GetParticleChangeForGamma();
00115 isInitialised = true;
00116 }
00117
00118
00119
00120 void G4LivermoreGammaConversionModel::ReadData(size_t Z, const char* path)
00121 {
00122 if (verboseLevel > 1)
00123 {
00124 G4cout << "Calling ReadData() of G4LivermoreGammaConversionModel"
00125 << G4endl;
00126 }
00127
00128 if(data[Z]) { return; }
00129
00130 const char* datadir = path;
00131
00132 if(!datadir)
00133 {
00134 datadir = getenv("G4LEDATA");
00135 if(!datadir)
00136 {
00137 G4Exception("G4LivermoreGammaConversionModel::ReadData()",
00138 "em0006",FatalException,
00139 "Environment variable G4LEDATA not defined");
00140 return;
00141 }
00142 }
00143
00144
00145
00146 data[Z] = new G4LPhysicsFreeVector();
00147
00148
00149 data[Z] ->SetSpline(true);
00150
00151
00152 std::ostringstream ost;
00153 ost << datadir << "/livermore/pair/pp-cs-" << Z <<".dat";
00154 std::ifstream fin(ost.str().c_str());
00155
00156 if( !fin.is_open())
00157 {
00158 G4ExceptionDescription ed;
00159 ed << "G4LivermoreGammaConversionModel data file <" << ost.str().c_str()
00160 << "> is not opened!" << G4endl;
00161 G4Exception("G4LivermoreGammaConversionModel::ReadData()",
00162 "em0003",FatalException,
00163 ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
00164 return;
00165 }
00166
00167 else
00168 {
00169
00170 if(verboseLevel > 3) { G4cout << "File " << ost.str()
00171 << " is opened by G4LivermoreGammaConversionModel" << G4endl;}
00172
00173 data[Z]->Retrieve(fin, true);
00174 }
00175
00176
00177 }
00178
00179
00180
00181 G4double
00182 G4LivermoreGammaConversionModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
00183 G4double GammaEnergy,
00184 G4double Z, G4double,
00185 G4double, G4double)
00186 {
00187 if (verboseLevel > 1)
00188 {
00189 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreGammaConversionModel"
00190 << G4endl;
00191 }
00192
00193 if (GammaEnergy < lowEnergyLimit) { return 0.0; }
00194
00195 G4double xs = 0.0;
00196
00197 G4int intZ=G4int(Z);
00198
00199 if(intZ < 1 || intZ > maxZ) { return xs; }
00200
00201 G4LPhysicsFreeVector* pv = data[intZ];
00202
00203
00204 if(!pv)
00205 {
00206 char* path = getenv("G4LEDATA");
00207 ReadData(intZ, path);
00208 pv = data[intZ];
00209 if(!pv) { return xs; }
00210 }
00211
00212 xs = pv->Value(GammaEnergy);
00213
00214 if(verboseLevel > 0)
00215 {
00216 G4int n = pv->GetVectorLength() - 1;
00217 G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)=" << GammaEnergy/MeV << G4endl;
00218 G4cout << " cs (Geant4 internal unit)=" << xs << G4endl;
00219 G4cout << " -> first cs value in EADL data file (iu) =" << (*pv)[0] << G4endl;
00220 G4cout << " -> last cs value in EADL data file (iu) =" << (*pv)[n] << G4endl;
00221 G4cout << "*********************************************************" << G4endl;
00222 }
00223
00224 return xs;
00225
00226 }
00227
00228
00229
00230 void G4LivermoreGammaConversionModel::SampleSecondaries(
00231 std::vector<G4DynamicParticle*>* fvect,
00232 const G4MaterialCutsCouple* couple,
00233 const G4DynamicParticle* aDynamicGamma,
00234 G4double, G4double)
00235 {
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247 if (verboseLevel > 1)
00248 G4cout << "Calling SampleSecondaries() of G4LivermoreGammaConversionModel" << G4endl;
00249
00250 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
00251 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
00252
00253 G4double epsilon ;
00254 G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
00255
00256
00257 if (photonEnergy < smallEnergy )
00258 {
00259 epsilon = epsilon0Local + (0.5 - epsilon0Local) * G4UniformRand();
00260 }
00261 else
00262 {
00263
00264
00265 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
00266 const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
00267
00268 if (element == 0)
00269 {
00270 G4cout << "G4LivermoreGammaConversionModel::SampleSecondaries - element = 0"
00271 << G4endl;
00272 return;
00273 }
00274 G4IonisParamElm* ionisation = element->GetIonisation();
00275 if (ionisation == 0)
00276 {
00277 G4cout << "G4LivermoreGammaConversionModel::SampleSecondaries - ionisation = 0"
00278 << G4endl;
00279 return;
00280 }
00281
00282
00283 G4double fZ = 8. * (ionisation->GetlogZ3());
00284 if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb());
00285
00286
00287 G4double screenFactor = 136. * epsilon0Local / (element->GetIonisation()->GetZ3()) ;
00288 G4double screenMax = std::exp ((42.24 - fZ)/8.368) - 0.952 ;
00289 G4double screenMin = std::min(4.*screenFactor,screenMax) ;
00290
00291
00292 G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
00293 G4double epsilonMin = std::max(epsilon0Local,epsilon1);
00294 G4double epsilonRange = 0.5 - epsilonMin ;
00295
00296
00297 G4double screen;
00298 G4double gReject ;
00299
00300 G4double f10 = ScreenFunction1(screenMin) - fZ;
00301 G4double f20 = ScreenFunction2(screenMin) - fZ;
00302 G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
00303 G4double normF2 = std::max(1.5 * f20,0.);
00304
00305 do
00306 {
00307 if (normF1 / (normF1 + normF2) > G4UniformRand() )
00308 {
00309 epsilon = 0.5 - epsilonRange * std::pow(G4UniformRand(), 0.333333) ;
00310 screen = screenFactor / (epsilon * (1. - epsilon));
00311 gReject = (ScreenFunction1(screen) - fZ) / f10 ;
00312 }
00313 else
00314 {
00315 epsilon = epsilonMin + epsilonRange * G4UniformRand();
00316 screen = screenFactor / (epsilon * (1 - epsilon));
00317 gReject = (ScreenFunction2(screen) - fZ) / f20 ;
00318 }
00319 } while ( gReject < G4UniformRand() );
00320
00321 }
00322
00323
00324
00325 G4double electronTotEnergy;
00326 G4double positronTotEnergy;
00327
00328 if (G4UniformRand() > 0.5)
00329 {
00330 electronTotEnergy = (1. - epsilon) * photonEnergy;
00331 positronTotEnergy = epsilon * photonEnergy;
00332 }
00333 else
00334 {
00335 positronTotEnergy = (1. - epsilon) * photonEnergy;
00336 electronTotEnergy = epsilon * photonEnergy;
00337 }
00338
00339
00340
00341
00342
00343 G4double u;
00344 const G4double a1 = 0.625;
00345 G4double a2 = 3. * a1;
00346
00347
00348
00349 if (0.25 > G4UniformRand())
00350 {
00351 u = - std::log(G4UniformRand() * G4UniformRand()) / a1 ;
00352 }
00353 else
00354 {
00355 u = - std::log(G4UniformRand() * G4UniformRand()) / a2 ;
00356 }
00357
00358 G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
00359 G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
00360 G4double phi = twopi * G4UniformRand();
00361
00362 G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
00363 G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
00364
00365
00366
00367
00368
00369
00370 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
00371
00372 G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
00373 electronDirection.rotateUz(photonDirection);
00374
00375 G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
00376 electronDirection,
00377 electronKineEnergy);
00378
00379
00380 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
00381
00382 G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
00383 positronDirection.rotateUz(photonDirection);
00384
00385
00386 G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
00387 positronDirection,
00388 positronKineEnergy);
00389
00390 fvect->push_back(particle1);
00391 fvect->push_back(particle2);
00392
00393
00394 fParticleChange->SetProposedKineticEnergy(0.);
00395 fParticleChange->ProposeTrackStatus(fStopAndKill);
00396
00397 }
00398
00399
00400
00401 G4double
00402 G4LivermoreGammaConversionModel::ScreenFunction1(G4double screenVariable)
00403 {
00404
00405
00406 G4double value;
00407
00408 if (screenVariable > 1.)
00409 value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
00410 else
00411 value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
00412
00413 return value;
00414 }
00415
00416
00417
00418 G4double
00419 G4LivermoreGammaConversionModel::ScreenFunction2(G4double screenVariable)
00420 {
00421
00422
00423 G4double value;
00424
00425 if (screenVariable > 1.)
00426 value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
00427 else
00428 value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
00429
00430 return value;
00431 }
00432