52 lowEnergyLimit = 2.0*electron_mass_c2;
62 G4cout <<
"G4LivermoreGammaConversionModel is constructed " <<
G4endl;
71 for(
G4int i=0; i<maxZ; ++i) {
88 G4cout <<
"Calling Initialise() of G4LivermoreGammaConversionModel."
91 << LowEnergyLimit() /
MeV <<
" MeV - "
92 << HighEnergyLimit() /
GeV <<
" GeV"
101 InitialiseElementSelectors(particle, cuts);
105 char* path = getenv(
"G4LEDATA");
112 for(
G4int i=0; i<numOfCouples; ++i)
119 for (
G4int j=0; j<nelm; ++j)
121 G4int Z = (
G4int)(*theElementVector)[j]->GetZ();
123 else if(Z > maxZ) { Z = maxZ; }
124 if(!data[Z]) { ReadData(Z, path); }
128 if(isInitialised) {
return; }
129 fParticleChange = GetParticleChangeForGamma();
130 isInitialised =
true;
148 return lowEnergyLimit;
155 if (verboseLevel > 1)
157 G4cout <<
"Calling ReadData() of G4LivermoreGammaConversionModel"
161 if(data[Z]) {
return; }
163 const char* datadir = path;
167 datadir = getenv(
"G4LEDATA");
170 G4Exception(
"G4LivermoreGammaConversionModel::ReadData()",
172 "Environment variable G4LEDATA not defined");
183 std::ostringstream ost;
184 ost << datadir <<
"/livermore/pair/pp-cs-" << Z <<
".dat";
185 std::ifstream fin(ost.str().c_str());
190 ed <<
"G4LivermoreGammaConversionModel data file <" << ost.str().c_str()
191 <<
"> is not opened!" <<
G4endl;
192 G4Exception(
"G4LivermoreGammaConversionModel::ReadData()",
194 ed,
"G4LEDATA version should be G4EMLOW6.27 or later.");
201 if(verboseLevel > 3) {
G4cout <<
"File " << ost.str()
202 <<
" is opened by G4LivermoreGammaConversionModel" <<
G4endl;}
204 data[Z]->Retrieve(fin,
true);
208 data[Z] ->SetSpline(
true);
220 if (verboseLevel > 1)
222 G4cout <<
"Calling ComputeCrossSectionPerAtom() of G4LivermoreGammaConversionModel"
226 if (GammaEnergy < lowEnergyLimit) {
return 0.0; }
232 if(intZ < 1 || intZ > maxZ) {
return xs; }
240 InitialiseForElement(0, intZ);
242 if(!pv) {
return xs; }
245 xs = pv->
Value(GammaEnergy);
250 G4cout <<
"****** DEBUG: tcs value for Z=" << Z <<
" at energy (MeV)="
252 G4cout <<
" cs (Geant4 internal unit)=" << xs <<
G4endl;
253 G4cout <<
" -> first cs value in EADL data file (iu) =" << (*pv)[0] <<
G4endl;
254 G4cout <<
" -> last cs value in EADL data file (iu) =" << (*pv)[
n] <<
G4endl;
255 G4cout <<
"*********************************************************" <<
G4endl;
265 std::vector<G4DynamicParticle*>* fvect,
281 if (verboseLevel > 1) {
282 G4cout <<
"Calling SampleSecondaries() of G4LivermoreGammaConversionModel"
290 G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
293 if (photonEnergy < smallEnergy )
295 epsilon = epsilon0Local + (0.5 - epsilon0Local) *
G4UniformRand();
302 const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
306 G4cout <<
"G4LivermoreGammaConversionModel::SampleSecondaries - element = 0"
313 G4cout <<
"G4LivermoreGammaConversionModel::SampleSecondaries - ionisation = 0"
320 if (photonEnergy > 50. *
MeV) fZ += 8. * (element->
GetfCoulomb());
328 G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
330 G4double epsilonRange = 0.5 - epsilonMin ;
336 G4double f10 = ScreenFunction1(screenMin) - fZ;
337 G4double f20 = ScreenFunction2(screenMin) - fZ;
345 epsilon = 0.5 - epsilonRange * std::pow(
G4UniformRand(), 0.333333) ;
346 screen = screenFactor / (epsilon * (1. - epsilon));
347 gReject = (ScreenFunction1(screen) - fZ) / f10 ;
352 screen = screenFactor / (epsilon * (1 - epsilon));
353 gReject = (ScreenFunction2(screen) - fZ) / f20 ;
366 electronTotEnergy = (1. - epsilon) * photonEnergy;
367 positronTotEnergy = epsilon * photonEnergy;
371 positronTotEnergy = (1. - epsilon) * photonEnergy;
372 electronTotEnergy = epsilon * photonEnergy;
394 G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
395 G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
398 G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
399 G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
406 G4double electronKineEnergy =
std::max(0.,electronTotEnergy - electron_mass_c2) ;
409 electronDirection.rotateUz(photonDirection);
416 G4double positronKineEnergy =
std::max(0.,positronTotEnergy - electron_mass_c2) ;
419 positronDirection.rotateUz(photonDirection);
426 fvect->push_back(particle1);
427 fvect->push_back(particle2);
430 fParticleChange->SetProposedKineticEnergy(0.);
444 if (screenVariable > 1.)
445 value = 42.24 - 8.368 *
G4Log(screenVariable + 0.952);
447 value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
461 if (screenVariable > 1.)
462 value = 42.24 - 8.368 *
G4Log(screenVariable + 0.952);
464 value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
478 G4AutoLock l(&LivermoreGammaConversionModelMutex);
481 if(!data[Z]) { ReadData(Z); }
virtual G4double MinPrimaryEnergy(const G4Material *, const G4ParticleDefinition *, G4double)
G4double ScreenFunction2(G4double screenVariable)
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
std::vector< G4Element * > G4ElementVector
std::ostringstream G4ExceptionDescription
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
void ReadData(size_t Z, const char *path=0)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
G4double GetfCoulomb() const
G4ParticleDefinition * GetDefinition() const
size_t GetVectorLength() const
const G4ElementVector * GetElementVector() const
G4double ScreenFunction1(G4double screenVariable)
#define G4MUTEX_INITIALIZER
static G4LPhysicsFreeVector * data[100]
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4GLOB_DLL std::ostream G4cout
size_t GetTableSize() const
G4LivermoreGammaConversionModel(const G4ParticleDefinition *p=0, const G4String &nam="LivermoreConversion")
const G4ThreeVector & GetMomentumDirection() const
virtual ~G4LivermoreGammaConversionModel()
std::vector< G4EmElementSelector * > * GetElementSelectors()
G4double Value(G4double theEnergy, size_t &lastidx) const
G4double GetlogZ3() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4double G4Log(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
static G4ProductionCutsTable * GetProductionCutsTable()
static G4Positron * Positron()
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4IonisParamElm * GetIonisation() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static G4Electron * Electron()
size_t GetNumberOfElements() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4ThreeVector G4ParticleMomentum
const G4Material * GetMaterial() const