63 theLorentzTables1(0),theLorentzTables2(0)
88 std::map<G4double,G4PhysicsTable*>::iterator j;
125 char* path = getenv(
"G4LEDATA");
129 "G4PenelopeBremsstrahlungAngular - G4LEDATA environment variable not set!";
130 G4Exception(
"G4PenelopeBremsstrahlungAngular::ReadDataFile()",
135 G4String pathFile = pathString +
"/penelope/bremsstrahlung/pdbrang.p08";
136 std::ifstream file(pathFile);
140 G4String excep =
"G4PenelopeBremsstrahlungAngular - data file " + pathFile +
" not found!";
141 G4Exception(
"G4PenelopeBremsstrahlungAngular::ReadDataFile()",
154 file >> iz1 >> ie1 >> ik1 >> zr >> er >> kr >> a1 >>
a2;
156 if ((iz1-1 == i) && (ik1-1 == k) && (ie1-1 == j))
164 ed <<
"Corrupted data file " << pathFile <<
"?" <<
G4endl;
165 G4Exception(
"G4PenelopeBremsstrahlungAngular::ReadDataFile()",
191 G4Exception(
"G4PenelopeBremsstrahlungAngular::PrepareInterpolationTables()",
202 const G4int reducedEnergyGrid=21;
225 QQ1vector->
PutValue(k,pZ[k],std::log(
QQ1[k][i][j]));
233 Q2[i][j]=QQ2vector->
Value(Zmat);
243 for(i=0;i<reducedEnergyGrid;i++)
248 betas[i]=std::sqrt(pE[i]*(pE[i]+2*electron_mass_c2))/(pE[i]+electron_mass_c2);
254 Q1[i][j]=Q1[i][j]/Zmat;
265 Q1vector->
PutValue(j,pK[j],std::log(Q1[i][j]));
266 Q2vector->
PutValue(j,pK[j],Q2[i][j]);
269 for (j=0;j<reducedEnergyGrid;j++)
271 Q1E[i][j]=Q1vector->
Value(ppK[j]);
272 Q2E[i][j]=Q2vector->
Value(ppK[j]);
288 for (j=0;j<reducedEnergyGrid;j++)
296 for (j=0;j<reducedEnergyGrid;j++)
302 thevec->
PutValue(i,betas[i],Q1E[i][j]);
303 thevec2->
PutValue(i,betas[i],Q2E[i][j]);
317 ed <<
"Unable to create tables of Lorentz coefficients for " <<
G4endl;
318 ed <<
"<Z>= " << Zmat <<
" in G4PenelopeBremsstrahlungAngular" <<
G4endl;
321 G4Exception(
"G4PenelopeBremsstrahlungAngular::PrepareInterpolationTables()",
337 G4Exception(
"G4PenelopeBremsstrahlungAngular::SampleDirection()",
347 G4Exception(
"G4PenelopeBremsstrahlungAngular::SampleDirection()",
357 G4Exception(
"G4PenelopeBremsstrahlungAngular::SampleDirection()",
358 "em2040",
FatalException,
"Material not found in the effectiveZ table");
364 G4cout <<
"Effective <Z> for material : " << material->
GetName() <<
370 G4double beta = std::sqrt(ePrimary*(ePrimary+2*electron_mass_c2))/
371 (ePrimary+electron_mass_c2);
377 if (ePrimary > 500*
keV)
383 cdt = -1.0*std::pow(-cdt,1./3.);
385 cdt = std::pow(cdt,1./3.);
387 cdt = (cdt+beta)/(1.0+beta*cdt);
389 sinTheta = std::sqrt(1. - cdt*cdt);
392 sinTheta* std::sin(phi),
403 ed <<
"Unable to retrieve Lorentz tables for Z= " << Zmat <<
G4endl;
404 G4Exception(
"G4PenelopeBremsstrahlungAngular::SampleDirection()",
421 P10 = v1->
Value(beta);
452 cdt = (cdt+betap)/(1.0+betap*cdt);
455 sinTheta = std::sqrt(1. - cdt*cdt);
458 sinTheta* std::sin(phi),
473 G4cout <<
"WARNING: G4PenelopeBremsstrahlungAngular() does NOT support PolarAngle()" <<
G4endl;
474 G4cout <<
"Please use the alternative interface SampleDirection()" <<
G4endl;
475 G4Exception(
"G4PenelopeBremsstrahlungAngular::PolarAngle()",
492 std::vector<G4double> *StechiometricFactors =
new std::vector<G4double>;
496 for (
G4int i=0;i<nElements;i++)
498 G4double fraction = fractionVector[i];
499 G4double atomicWeigth = (*elementVector)[i]->GetA()/(
g/
mole);
500 StechiometricFactors->push_back(fraction/atomicWeigth);
503 G4double MaxStechiometricFactor = 0.;
504 for (
G4int i=0;i<nElements;i++)
506 if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
507 MaxStechiometricFactor = (*StechiometricFactors)[i];
510 for (
G4int i=0;i<nElements;i++)
511 (*StechiometricFactors)[i] /= MaxStechiometricFactor;
515 for (
G4int i=0;i<nElements;i++)
517 G4double Z = (*elementVector)[i]->GetZ();
518 sumz2 += (*StechiometricFactors)[i]*Z*Z;
519 sums += (*StechiometricFactors)[i];
521 delete StechiometricFactors;
523 G4double ZBR = std::sqrt(sumz2/sums);
std::map< G4double, G4PhysicsTable * > * theLorentzTables2
static const G4double * P1[nN]
std::vector< G4Element * > G4ElementVector
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
std::ostringstream G4ExceptionDescription
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
const G4String & GetName() const
void push_back(G4PhysicsVector *)
static const G4double P20[nE]
static const G4double P10[nE]
static const G4double P11[nE]
void PrepareTables(const G4Material *material, G4bool isMaster)
Reserved for Master Model.
G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double out_energy, G4int Z, const G4Material *mat=0)
Samples the direction of the outgoing photon (in global coordinates).
const G4ElementVector * GetElementVector() const
std::map< G4double, G4PhysicsTable * > * theLorentzTables1
static const G4int NumberofKPoints
G4GLOB_DLL std::ostream G4cout
std::ostream & tab(std::ostream &)
const G4ThreeVector & GetMomentumDirection() const
G4double QQ1[NumberofZPoints][NumberofEPoints][NumberofKPoints]
static const double twopi
G4double Value(G4double theEnergy, size_t &lastidx) const
void Initialize()
Reserved for Master Model The Initialize() method forces the cleaning of tables.
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double CalculateEffectiveZ(const G4Material *material)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static const G4double * P2[nN]
G4double QQ2[NumberofZPoints][NumberofEPoints][NumberofKPoints]
G4PenelopeBremsstrahlungAngular()
size_t GetNumberOfElements() const
static const G4int NumberofEPoints
G4ThreeVector fLocalDirection
G4double PolarAngle(const G4double initial_energy, const G4double final_energy, const G4int Z)
Old interface, backwards compatibility. Will not work in this case it will produce a G4Exception()...
std::map< const G4Material *, G4double > * theEffectiveZSq
const G4double * GetFractionVector() const
static const G4int NumberofZPoints
~G4PenelopeBremsstrahlungAngular()
static const G4double P21[nE]