43G4int G4BoldyshevTripletModel::maxZ = 99;
49 fParticleChange =
nullptr;
51 lowEnergyLimit = 4.0*electron_mass_c2;
52 momentumThreshold_c = energyThreshold = xb = xn = lowEnergyLimit;
62 G4cout <<
"G4BoldyshevTripletModel is constructed " <<
G4endl;
71 for(
G4int i=0; i<maxZ; ++i) {
87 G4cout <<
"Calling Initialise() of G4BoldyshevTripletModel."
95 energyThreshold = 1.1*electron_mass_c2;
96 momentumThreshold_c = std::sqrt(energyThreshold * energyThreshold
97 - electron_mass_c2*electron_mass_c2);
98 G4double momentumThreshold_N = momentumThreshold_c/electron_mass_c2;
100 std::sqrt(momentumThreshold_N*momentumThreshold_N + 1.0));
105 G4double J1 = 0.5*(t*cosht/sinht - logsinht);
106 G4double J2 = (-2./3.)*logsinht + t*cosht/sinht
107 + (sinht - t*cosht*cosht*cosht)/(3.*sinht*sinht*sinht);
115 char* path = std::getenv(
"G4LEDATA");
122 for(
G4int i=0; i<numOfCouples; ++i)
129 for (
G4int j=0; j<nelm; ++j)
131 G4int Z = std::min((*theElementVector)[j]->GetZasInt(), maxZ);
132 if(!data[Z]) { ReadData(Z, path); }
136 if(!fParticleChange) {
148 return lowEnergyLimit;
153void G4BoldyshevTripletModel::ReadData(
size_t Z,
const char* path)
155 if (verboseLevel > 1)
157 G4cout <<
"Calling ReadData() of G4BoldyshevTripletModel"
161 if(data[Z]) {
return; }
163 const char* datadir = path;
167 datadir = std::getenv(
"G4LEDATA");
172 "Environment variable G4LEDATA not defined");
178 std::ostringstream ost;
179 ost << datadir <<
"/livermore/tripdata/pp-trip-cs-" << Z <<
".dat";
180 std::ifstream fin(ost.str().c_str());
185 ed <<
"G4BoldyshevTripletModel data file <" << ost.str().c_str()
186 <<
"> is not opened!" <<
G4endl;
189 ed,
"G4LEDATA version should be G4EMLOW6.27 or later.");
196 if(verboseLevel > 3) {
G4cout <<
"File " << ost.str()
197 <<
" is opened by G4BoldyshevTripletModel" <<
G4endl;}
212 if (verboseLevel > 1)
214 G4cout <<
"Calling ComputeCrossSectionPerAtom() of G4BoldyshevTripletModel"
218 if (GammaEnergy < lowEnergyLimit) {
return 0.0; }
230 if(!pv) {
return xs; }
233 xs = pv->
Value(GammaEnergy);
237 G4cout <<
"*** Triplet conversion xs for Z=" << Z <<
" at energy E(MeV)="
238 << GammaEnergy/MeV <<
" cs=" << xs/millibarn <<
" mb" <<
G4endl;
246 std::vector<G4DynamicParticle*>* fvect,
254 if (verboseLevel > 1) {
255 G4cout <<
"Calling SampleSecondaries() of G4BoldyshevTripletModel"
268 static const G4double costlim = std::cos(4.47*CLHEP::pi/180.);
270 G4double loga, f1_re, greject, cost;
271 G4double cosThetaMax = (energyThreshold - electron_mass_c2
272 + electron_mass_c2*(energyThreshold + electron_mass_c2)/photonEnergy )
273 /momentumThreshold_c;
274 if (cosThetaMax > 1.) {
283 cost =
G4Exp(logcostm*rndmEngine->
flat());
285 G4double bre = (1.-5.*cost*cost)/(2.*cost);
286 loga =
G4Log((1.+ cost)/(1.- cost));
287 f1_re = 1. - bre*loga;
288 greject = (cost < costlim) ? are*f1_re : 1.0;
291 }
while(greject < rndmEngine->flat());
294 G4double sint2 = (1. - cost)*(1. + cost);
295 G4double fp = 1. - sint2*loga/(2.*cost) ;
299 phi_re = twopi*rndmEngine->
flat();
300 rt = (1. - std::cos(2.*phi_re)*fp/f1_re)/twopi;
303 }
while(rt < rndmEngine->flat());
306 G4double S = electron_mass_c2*(2.* photonEnergy + electron_mass_c2);
307 G4double P2 =
S - electron_mass_c2*electron_mass_c2;
309 G4double D2 = 4.*
S * electron_mass_c2*electron_mass_c2 + P2*P2*sint2;
310 G4double ener_re = electron_mass_c2 * (
S + electron_mass_c2*electron_mass_c2)/sqrt(D2);
312 if(ener_re >= energyThreshold)
314 G4double electronRKineEnergy = ener_re - electron_mass_c2;
316 G4ThreeVector electronRDirection (sint*std::cos(phi_re), sint*std::sin(phi_re), cost);
317 electronRDirection.
rotateUz(photonDirection);
320 electronRKineEnergy);
335 G4double a = std::sqrt(16./xb - 3. - 36.*re*xn + 36.*re*re*xn*xn + 6.*xb*re*xn);
337 epsilon = c1/(2.*xb) + (xb - 4.)/(2.*c1) + 0.5;
339 G4double photonEnergy1 = photonEnergy - ener_re ;
341 G4double positronTotEnergy = std::max(
epsilon*photonEnergy1, electron_mass_c2);
342 G4double electronTotEnergy = std::max(photonEnergy1 - positronTotEnergy, electron_mass_c2);
345 static const G4double a2 = 0.5333333333;
347 G4double u = (0.25 > rndmEngine->
flat()) ? uu*a1 : uu*a2;
349 G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
350 G4double sinte = std::sin(thetaEle);
351 G4double coste = std::cos(thetaEle);
353 G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
354 G4double sintp = std::sin(thetaPos);
355 G4double costp = std::cos(thetaPos);
365 G4double electronKineEnergy = electronTotEnergy - electron_mass_c2;
367 G4ThreeVector electronDirection (sinte*cosp, sinte*sinp, coste);
368 electronDirection.
rotateUz(photonDirection);
374 G4double positronKineEnergy = positronTotEnergy - electron_mass_c2;
376 G4ThreeVector positronDirection (-sintp*cosp, -sintp*sinp, costp);
377 positronDirection.
rotateUz(photonDirection);
381 positronDirection, positronKineEnergy);
384 fvect->push_back(particle1);
385 fvect->push_back(particle2);
387 if(particle3) { fvect->push_back(particle3); }
405 if(!data[Z]) { ReadData(Z); }
double epsilon(double density, double temperature)
std::vector< G4Element * > G4ElementVector
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double G4Log(G4double x)
#define G4MUTEX_INITIALIZER
G4GLOB_DLL std::ostream G4cout
Hep3Vector & rotateUz(const Hep3Vector &)
virtual G4double MinPrimaryEnergy(const G4Material *, const G4ParticleDefinition *, G4double)
virtual ~G4BoldyshevTripletModel()
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cut=0., G4double emax=DBL_MAX)
G4BoldyshevTripletModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="BoldyshevTripletConversion")
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
const G4Material * GetMaterial() const
const G4ElementVector * GetElementVector() const
size_t GetNumberOfElements() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4double Value(G4double theEnergy, std::size_t &lastidx) const
virtual G4bool Retrieve(std::ifstream &fIn, G4bool ascii=false)
static G4Positron * Positron()
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
std::size_t GetTableSize() const
static G4ProductionCutsTable * GetProductionCutsTable()
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4double LowEnergyLimit() const
G4double HighEnergyLimit() const
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)