43 G4int G4LivermoreGammaConversionModel::maxZ = 99;
62 G4cout <<
"G4LivermoreGammaConversionModel is constructed " <<
G4endl;
79 G4cout <<
"Calling Initialise() of G4LivermoreGammaConversionModel."
82 << LowEnergyLimit() /
MeV <<
" MeV - "
83 << HighEnergyLimit() /
GeV <<
" GeV"
92 InitialiseElementSelectors(particle, cuts);
96 char* path = getenv(
"G4LEDATA");
103 for(
G4int i=0; i<numOfCouples; ++i)
110 for (
G4int j=0; j<nelm; ++j)
114 else if(Z > maxZ) { Z = maxZ; }
115 if(!
data[Z]) { ReadData(Z, path); }
122 if(isInitialised) {
return; }
123 fParticleChange = GetParticleChangeForGamma();
124 isInitialised =
true;
142 return lowEnergyLimit;
147 void G4LivermoreGammaConversionModel::ReadData(
size_t Z,
const char* path)
149 if (verboseLevel > 1)
151 G4cout <<
"Calling ReadData() of G4LivermoreGammaConversionModel"
155 if(
data[Z]) {
return; }
157 const char* datadir = path;
161 datadir = getenv(
"G4LEDATA");
164 G4Exception(
"G4LivermoreGammaConversionModel::ReadData()",
166 "Environment variable G4LEDATA not defined");
177 std::ostringstream ost;
178 ost << datadir <<
"/livermore/pair/pp-cs-" << Z <<
".dat";
179 std::ifstream
fin(ost.str().c_str());
184 ed <<
"G4LivermoreGammaConversionModel data file <" << ost.str().c_str()
185 <<
"> is not opened!" <<
G4endl;
186 G4Exception(
"G4LivermoreGammaConversionModel::ReadData()",
188 ed,
"G4LEDATA version should be G4EMLOW6.27 or later.");
195 if(verboseLevel > 3) {
G4cout <<
"File " << ost.str()
196 <<
" is opened by G4LivermoreGammaConversionModel" <<
G4endl;}
202 data[
Z] ->SetSpline(
true);
214 if (verboseLevel > 1)
216 G4cout <<
"Calling ComputeCrossSectionPerAtom() of G4LivermoreGammaConversionModel"
220 if (GammaEnergy < lowEnergyLimit) {
return 0.0; }
226 if(intZ < 1 || intZ > maxZ) {
return xs; }
234 InitialiseForElement(0, intZ);
236 if(!pv) {
return xs; }
239 xs = pv->
Value(GammaEnergy);
244 G4cout <<
"****** DEBUG: tcs value for Z=" << Z <<
" at energy (MeV)="
246 G4cout <<
" cs (Geant4 internal unit)=" << xs <<
G4endl;
247 G4cout <<
" -> first cs value in EADL data file (iu) =" << (*pv)[0] <<
G4endl;
248 G4cout <<
" -> last cs value in EADL data file (iu) =" << (*pv)[
n] <<
G4endl;
249 G4cout <<
"*********************************************************" <<
G4endl;
259 std::vector<G4DynamicParticle*>* fvect,
275 if (verboseLevel > 1) {
276 G4cout <<
"Calling SampleSecondaries() of G4LivermoreGammaConversionModel"
287 if (photonEnergy < smallEnergy )
289 epsilon = epsilon0Local + (0.5 - epsilon0Local) *
G4UniformRand();
296 const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
300 G4cout <<
"G4LivermoreGammaConversionModel::SampleSecondaries - element = 0"
307 G4cout <<
"G4LivermoreGammaConversionModel::SampleSecondaries - ionisation = 0"
314 if (photonEnergy > 50. *
MeV) fZ += 8. * (element->
GetfCoulomb());
322 G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
324 G4double epsilonRange = 0.5 - epsilonMin ;
330 G4double f10 = ScreenFunction1(screenMin) - fZ;
331 G4double f20 = ScreenFunction2(screenMin) - fZ;
339 epsilon = 0.5 - epsilonRange * std::pow(
G4UniformRand(), 0.333333) ;
340 screen = screenFactor / (epsilon * (1. - epsilon));
341 gReject = (ScreenFunction1(screen) - fZ) / f10 ;
346 screen = screenFactor / (epsilon * (1 - epsilon));
347 gReject = (ScreenFunction2(screen) - fZ) / f20 ;
360 electronTotEnergy = (1. - epsilon) * photonEnergy;
361 positronTotEnergy = epsilon * photonEnergy;
365 positronTotEnergy = (1. - epsilon) * photonEnergy;
366 electronTotEnergy = epsilon * photonEnergy;
392 G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
393 G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
403 electronDirection.
rotateUz(photonDirection);
413 positronDirection.
rotateUz(photonDirection);
420 fvect->push_back(particle1);
421 fvect->push_back(particle2);
424 fParticleChange->SetProposedKineticEnergy(0.);
432 G4LivermoreGammaConversionModel::ScreenFunction1(
G4double screenVariable)
438 if (screenVariable > 1.)
439 value = 42.24 - 8.368 *
G4Log(screenVariable + 0.952);
441 value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
449 G4LivermoreGammaConversionModel::ScreenFunction2(
G4double screenVariable)
455 if (screenVariable > 1.)
456 value = 42.24 - 8.368 *
G4Log(screenVariable + 0.952);
458 value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
472 G4AutoLock l(&LivermoreGammaConversionModelMutex);
475 if(!
data[Z]) { ReadData(Z); }
virtual G4double MinPrimaryEnergy(const G4Material *, const G4ParticleDefinition *, G4double)
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
std::vector< G4Element * > G4ElementVector
std::ostringstream G4ExceptionDescription
G4double GetKineticEnergy() const
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
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
#define G4MUTEX_INITIALIZER
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
Hep3Vector & rotateUz(const Hep3Vector &)
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
const XML_Char int const XML_Char * value
static G4Electron * Electron()
size_t GetNumberOfElements() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
const XML_Char const XML_Char * data
const G4Material * GetMaterial() const