54 theReducedXSTable(0),theEffectiveZSq(0),theSamplingTable(0),
55 thePBcut(0),fVerbosity(verbosity)
59 {1.0e-12,0.025e0,0.05e0,0.075e0,0.1e0,0.15e0,0.2e0,0.25e0,
60 0.3e0,0.35e0,0.4e0,0.45e0,0.5e0,0.55e0,0.6e0,0.65e0,0.7e0,
61 0.75e0,0.8e0,0.85e0,0.9e0,0.925e0,0.95e0,0.97e0,0.99e0,
62 0.995e0,0.999e0,0.9995e0,0.9999e0,0.99995e0,0.99999e0,1.0e0};
64 for (
size_t ix=0;ix<
nBinsX;ix++)
67 for (
size_t i=0;i<
nBinsE;i++)
100 G4Exception(
"G4PenelopeBremsstrahlungFS::ClearTables()",
155 ed <<
"The container for the <Z^2> values is not initialized" <<
G4endl;
156 G4Exception(
"G4PenelopeBremsstrahlungFS::GetEffectiveZSquared()",
166 ed <<
"The value of <Z^2> is not properly set for material " <<
169 G4Exception(
"G4PenelopeBremsstrahlungFS::GetEffectiveZSquared()",
189 G4Exception(
"G4PenelopeBremsstrahlungFS::BuildScaledXSTable()",
194 G4cout <<
"Entering in G4PenelopeBremsstrahlungFS::BuildScaledXSTable for " <<
196 G4cout <<
"Threshold = " << cut/
keV <<
" keV, isMaster= " << isMaster <<
203 new std::map< std::pair<const G4Material*,G4double> ,
G4PhysicsTable*>;
220 std::vector<G4double> *StechiometricFactors =
new std::vector<G4double>;
224 for (
G4int i=0;i<nElements;i++)
226 G4double fraction = fractionVector[i];
227 G4double atomicWeigth = (*elementVector)[i]->GetA()/(
g/
mole);
228 StechiometricFactors->push_back(fraction/atomicWeigth);
231 G4double MaxStechiometricFactor = 0.;
232 for (
G4int i=0;i<nElements;i++)
234 if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
235 MaxStechiometricFactor = (*StechiometricFactors)[i];
238 for (
G4int i=0;i<nElements;i++)
239 (*StechiometricFactors)[i] /= MaxStechiometricFactor;
243 for (
G4int i=0;i<nElements;i++)
245 G4double Z = (*elementVector)[i]->GetZ();
246 sumz2 += (*StechiometricFactors)[i]*Z*
Z;
247 sums += (*StechiometricFactors)[i];
260 for (
G4int iel=0;iel<nElements;iel++)
262 G4double Z = (*elementVector)[iel]->GetZ();
264 G4double wgt = (*StechiometricFactors)[iel]*Z*Z/ZBR2;
274 ed <<
"Error in G4PenelopeBremsstrahlungFS::BuildScaledXSTable" <<
G4endl;
275 ed <<
"Unable to retrieve data for element " << iZ <<
G4endl;
276 G4Exception(
"G4PenelopeBremsstrahlungFS::BuildScaledXSTable()",
283 for (
size_t ie=0;ie<
nBinsE;ie++)
285 (*tempData)[ie] += wgt*(*atomData)[ie*(nBinsX+1)+nBinsX];
286 for (
size_t ix=0;ix<
nBinsX;ix++)
287 (*tempMatrix)[ie*nBinsX+ix] += wgt*(*atomData)[ie*(nBinsX+1)+ix];
295 for (
size_t ie=0;ie<
nBinsE;ie++)
299 for (
size_t ix=0;ix<
nBinsX;ix++)
300 tempData2[ix] = (*tempMatrix)[ie*nBinsX+ix];
305 G4double fnorm = (*tempData)[ie]/(rsum*fact);
306 G4double TST = 100.*std::fabs(fnorm-1.0);
310 ed <<
"G4PenelopeBremsstrahlungFS. Corrupted data files?" <<
G4endl;
311 G4cout <<
"TST= " << TST <<
"; fnorm = " << fnorm <<
G4endl;
315 G4Exception(
"G4PenelopeBremsstrahlungFS::BuildScaledXSTable()",
318 for (
size_t ix=0;ix<
nBinsX;ix++)
319 (*tempMatrix)[ie*nBinsX+ix] *= fnorm;
325 G4PhysicsTable* thePhysicsTable =
new G4PhysicsTable();
333 for (
size_t i=0;i<
nBinsX;i++)
334 thePhysicsTable->
push_back(
new G4PhysicsFreeVector(nBinsE+1));
336 for (
size_t ix=0;ix<
nBinsX;ix++)
338 G4PhysicsFreeVector* theVec =
339 (G4PhysicsFreeVector*) ((*thePhysicsTable)[ix]);
340 for (
size_t ie=0;ie<
nBinsE;ie++)
343 G4double aValue = (*tempMatrix)[ie*nBinsX+ix];
346 theVec->
PutValue(ie+1,logene,std::log(aValue));
352 G4double val1eV = (*theVec)[1]+derivative*(log1eV-theVec->
Energy(1));
356 std::pair<const G4Material*,G4double> theKey = std::make_pair(material,cut);
359 delete StechiometricFactors;
375 char* path = getenv(
"G4LEDATA");
378 G4String excep =
"G4PenelopeBremsstrahlungFS - G4LEDATA environment variable not set!";
379 G4Exception(
"G4PenelopeBremsstrahlungFS::ReadDataFile()",
386 std::ostringstream ost;
388 ost << path <<
"/penelope/bremsstrahlung/pdebr" << Z <<
".p08";
390 ost << path <<
"/penelope/bremsstrahlung/pdebr0" << Z <<
".p08";
391 std::ifstream
file(ost.str().c_str());
394 G4String excep =
"G4PenelopeBremsstrahlungFS - data file " +
395 G4String(ost.str()) +
" not found!";
396 G4Exception(
"G4PenelopeBremsstrahlungFS::ReadDataFile()",
408 ed <<
"Corrupted data file for Z=" << Z <<
G4endl;
409 G4Exception(
"G4PenelopeBremsstrahlungFS::ReadDataFile()",
416 for (
size_t ie=0;ie<
nBinsE;ie++)
423 for (
size_t ix=0;ix<
nBinsX;ix++)
426 (*theMatrix)[ie*(nBinsX+1)+ix] = myDouble*
millibarn;
429 (*theMatrix)[ie*(nBinsX+1)+nBinsX] = myDouble*
millibarn;
450 size_t size = nBinsX;
454 if (momOrder<-1 || size<2 || theXGrid[0]<0)
456 G4Exception(
"G4PenelopeBremsstrahlungFS::GetMomentumIntegral()",
460 for (
size_t i=1;i<size;i++)
462 if (theXGrid[i]<0 || theXGrid[i]<theXGrid[i-1])
465 ed <<
"Invalid call for bin " << i <<
G4endl;
466 G4Exception(
"G4PenelopeBremsstrahlungFS::GetMomentumIntegral()",
473 if (xup < theXGrid[0])
475 bool loopAgain =
true;
478 for (
size_t i=0;i<size-1;i++)
494 if (std::fabs(dx)>1
e-14*std::fabs(dy))
499 ds = a*std::log(xtc/x1)+b*(xtc-
x1);
500 else if (momOrder == 0)
501 ds = a*(xtc-
x1) + 0.5*b*(xtc*xtc-x1*x1);
503 ds = a*(std::pow(xtc,momOrder+1)-std::pow(x1,momOrder+1))/((
G4double) (momOrder + 1))
504 + b*(std::pow(xtc,momOrder+2)-std::pow(x1,momOrder+2))/((
G4double) (momOrder + 2));
507 ds = 0.5*(y1+
y2)*(xtc-x1)*std::pow(xtc,momOrder);
521 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
525 G4Exception(
"G4PenelopeBremsstrahlungFS::GetScaledXSTable()",
526 "em2013",
FatalException,
"Unable to retrieve the cross section table");
538 G4cout <<
"Entering in G4PenelopeBremsstrahlungFS::InitializeEnergySampling() for " <<
542 std::pair<const G4Material*,G4double> theKey = std::make_pair(material,cut);
550 for (
size_t i=0;i<
nBinsE;i++)
556 G4Exception(
"G4PenelopeBremsstrahlungFS::InitializeEnergySampling()",
557 "em2013",
FatalException,
"Unable to retrieve the cross section table");
560 for (
size_t ie=0;ie<
nBinsE;ie++)
567 for (
size_t ix=1;ix<
nBinsX;ix++)
591 for (
size_t ix=0;ix<
nBinsX;ix++)
594 tempData[ix] =
G4Exp((*vv)[ie+1]);
604 thePBcut->insert(std::make_pair(theKey,thePBvec));
613 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
617 ed <<
"Unable to retrieve the SamplingTable: " <<
620 G4Exception(
"G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
630 G4bool firstOrLastBin =
false;
635 firstOrLastBin =
true;
640 firstOrLastBin =
true;
673 G4cout <<
"Creating new instance of G4PhysicsFreeVector() on the worker" <<
G4endl;
681 for (
size_t iloop=0;iloop<
nBinsX;iloop++)
683 G4double val = (*theVec1)[iloop]+(((*theVec2)[iloop]-(*theVec1)[iloop]))*
690 for (
size_t iloop=0;iloop<
nBinsX;iloop++)
699 pbcut = (*(
thePBcut->find(theKey)->second))[eBin] +
700 ((*(
thePBcut->find(theKey)->second))[eBin+1]-(*(
thePBcut->find(theKey)->second))[eBin])*
707 G4int nIterations = 0;
715 if (pt < (*theTempVec)[0])
717 else if (pt > (*theTempVec)[
nBinsX-1])
722 if (delta < pt*1
e-10)
726 ed <<
"Found that (pt > (*theTempVec)[nBinsX-1]) with pt = " << pt <<
727 " , (*theTempVec)[nBinsX-1] = " << (*theTempVec)[
nBinsX-1] <<
" and delta = " <<
729 ed <<
"Possible symptom of problem with numerical precision" <<
G4endl;
730 G4Exception(
"G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
736 ed <<
"Crash at (pt > (*theTempVec)[nBinsX-1]) with pt = " << pt <<
737 " , (*theTempVec)[nBinsX-1]=" << (*theTempVec)[
nBinsX-1] <<
" and nBinsX = " <<
739 ed <<
"Material: " << mat->
GetName() <<
", energy = " << energy/
keV <<
" keV" <<
741 G4Exception(
"G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
752 if (pt > (*theTempVec)[k])
790 ed <<
"Warning in G4PenelopeBremsstrahlungFS::SampleX()" <<
G4endl;
791 ed <<
"Conflicting end-point values: w1=" << w1 <<
"; w2 = " << w2 <<
G4endl;
792 ed <<
"wbcut = " << wbcut <<
" energy= " << energy/
keV <<
" keV" <<
G4endl;
793 ed <<
"cut = " << cut/
keV <<
" keV" <<
G4endl;
794 G4Exception(
"G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
"em2015",
810 if (nIterations > 100)
812 }
while(eGamma < cut);
G4double Energy(size_t index) const
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const G4double * GetFractionVector() const
std::ostringstream G4ExceptionDescription
std::vector< ExP01TrackerHit * > a
G4double theXGrid[nBinsX]
G4Cache< G4PhysicsFreeVector * > fCache
Float_t y1[n_points_granero]
static constexpr double keV
void ReadDataFile(G4int Z)
Float_t x1[n_points_granero]
void InitializeEnergySampling(const G4Material *, G4double cut)
G4double GetEffectiveZSquared(const G4Material *mat) const
void BuildScaledXSTable(const G4Material *material, G4double cut, G4bool isMaster)
const G4ElementVector * GetElementVector() const
G4double theEGrid[nBinsE]
std::map< std::pair< const G4Material *, G4double >, G4PhysicsTable * > * theReducedXSTable
static constexpr double classic_electr_radius
std::map< std::pair< const G4Material *, G4double >, G4PhysicsTable * > * theSamplingTable
~G4PenelopeBremsstrahlungFS()
Float_t y2[n_points_geant4]
const G4String & GetName() const
static constexpr double g
G4PenelopeBremsstrahlungFS(G4int verbosity=0)
Only master models are supposed to create instances.
static constexpr double millibarn
const XML_Char int const XML_Char * value
static constexpr double electron_mass_c2
G4double GetMomentumIntegral(G4double *y, G4double up, G4int momOrder) const
std::map< std::pair< const G4Material *, G4double >, G4PhysicsFreeVector * > * thePBcut
double A(double temperature)
static constexpr double eV
static const size_t nBinsE
void ClearTables(G4bool isMaster=true)
Reserved for the master model: they build and handle tables.
void Put(const value_type &val) const
G4double SampleGammaEnergy(G4double energy, const G4Material *, const G4double cut) const
G4double G4ParticleHPJENDLHEData::G4double result
std::map< const G4Material *, G4double > * theEffectiveZSq
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::map< G4int, G4DataVector * > * theElementData
void push_back(G4PhysicsVector *)
std::vector< G4Element * > G4ElementVector
static constexpr double barn
G4GLOB_DLL std::ostream G4cout
const G4PhysicsTable * GetScaledXSTable(const G4Material *, const G4double cut) const
void PutValue(size_t index, G4double energy, G4double dataValue)
static constexpr double mole
static const G4double eps
static constexpr double fine_structure_const
Float_t x2[n_points_geant4]
double B(double temperature)
size_t GetNumberOfElements() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static const size_t nBinsX