53 theReducedXSTable(0),theEffectiveZSq(0),theSamplingTable(0),
54 thePBcut(0),fVerbosity(verbosity)
58 {1.0e-12,0.025e0,0.05e0,0.075e0,0.1e0,0.15e0,0.2e0,0.25e0,
59 0.3e0,0.35e0,0.4e0,0.45e0,0.5e0,0.55e0,0.6e0,0.65e0,0.7e0,
60 0.75e0,0.8e0,0.85e0,0.9e0,0.925e0,0.95e0,0.97e0,0.99e0,
61 0.995e0,0.999e0,0.9995e0,0.9999e0,0.99995e0,0.99999e0,1.0e0};
63 for (
size_t ix=0;ix<
nBinsX;ix++)
66 for (
size_t i=0;i<
nBinsE;i++)
99 G4Exception(
"G4PenelopeBremsstrahlungFS::ClearTables()",
154 ed <<
"The container for the <Z^2> values is not initialized" <<
G4endl;
155 G4Exception(
"G4PenelopeBremsstrahlungFS::GetEffectiveZSquared()",
165 ed <<
"The value of <Z^2> is not properly set for material " <<
168 G4Exception(
"G4PenelopeBremsstrahlungFS::GetEffectiveZSquared()",
188 G4Exception(
"G4PenelopeBremsstrahlungFS::BuildScaledXSTable()",
193 G4cout <<
"Entering in G4PenelopeBremsstrahlungFS::BuildScaledXSTable for " <<
195 G4cout <<
"Threshold = " << cut/
keV <<
" keV, isMaster= " << isMaster <<
202 new std::map< std::pair<const G4Material*,G4double> ,
G4PhysicsTable*>;
219 std::vector<G4double> *StechiometricFactors =
new std::vector<G4double>;
223 for (
G4int i=0;i<nElements;i++)
225 G4double fraction = fractionVector[i];
226 G4double atomicWeigth = (*elementVector)[i]->GetA()/(
g/
mole);
227 StechiometricFactors->push_back(fraction/atomicWeigth);
230 G4double MaxStechiometricFactor = 0.;
231 for (
G4int i=0;i<nElements;i++)
233 if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
234 MaxStechiometricFactor = (*StechiometricFactors)[i];
237 for (
G4int i=0;i<nElements;i++)
238 (*StechiometricFactors)[i] /= MaxStechiometricFactor;
242 for (
G4int i=0;i<nElements;i++)
244 G4double Z = (*elementVector)[i]->GetZ();
245 sumz2 += (*StechiometricFactors)[i]*Z*
Z;
246 sums += (*StechiometricFactors)[i];
259 for (
G4int iel=0;iel<nElements;iel++)
261 G4double Z = (*elementVector)[iel]->GetZ();
263 G4double wgt = (*StechiometricFactors)[iel]*Z*Z/ZBR2;
273 ed <<
"Error in G4PenelopeBremsstrahlungFS::BuildScaledXSTable" <<
G4endl;
274 ed <<
"Unable to retrieve data for element " << iZ <<
G4endl;
275 G4Exception(
"G4PenelopeBremsstrahlungFS::BuildScaledXSTable()",
282 for (
size_t ie=0;ie<
nBinsE;ie++)
284 (*tempData)[ie] += wgt*(*atomData)[ie*(nBinsX+1)+nBinsX];
285 for (
size_t ix=0;ix<
nBinsX;ix++)
286 (*tempMatrix)[ie*nBinsX+ix] += wgt*(*atomData)[ie*(nBinsX+1)+ix];
294 for (
size_t ie=0;ie<
nBinsE;ie++)
298 for (
size_t ix=0;ix<
nBinsX;ix++)
299 tempData2[ix] = (*tempMatrix)[ie*nBinsX+ix];
304 G4double fnorm = (*tempData)[ie]/(rsum*fact);
305 G4double TST = 100.*std::fabs(fnorm-1.0);
309 ed <<
"G4PenelopeBremsstrahlungFS. Corrupted data files?" <<
G4endl;
310 G4cout <<
"TST= " << TST <<
"; fnorm = " << fnorm <<
G4endl;
314 G4Exception(
"G4PenelopeBremsstrahlungFS::BuildScaledXSTable()",
317 for (
size_t ix=0;ix<
nBinsX;ix++)
318 (*tempMatrix)[ie*nBinsX+ix] *= fnorm;
324 G4PhysicsTable* thePhysicsTable =
new G4PhysicsTable();
332 for (
size_t i=0;i<
nBinsX;i++)
333 thePhysicsTable->
push_back(
new G4PhysicsFreeVector(nBinsE+1));
335 for (
size_t ix=0;ix<
nBinsX;ix++)
337 G4PhysicsFreeVector* theVec =
338 (G4PhysicsFreeVector*) ((*thePhysicsTable)[ix]);
339 for (
size_t ie=0;ie<
nBinsE;ie++)
342 G4double aValue = (*tempMatrix)[ie*nBinsX+ix];
345 theVec->
PutValue(ie+1,logene,std::log(aValue));
351 G4double val1eV = (*theVec)[1]+derivative*(log1eV-theVec->
Energy(1));
355 std::pair<const G4Material*,G4double> theKey = std::make_pair(material,cut);
358 delete StechiometricFactors;
374 char* path = std::getenv(
"G4LEDATA");
377 G4String excep =
"G4PenelopeBremsstrahlungFS - G4LEDATA environment variable not set!";
378 G4Exception(
"G4PenelopeBremsstrahlungFS::ReadDataFile()",
385 std::ostringstream ost;
387 ost << path <<
"/penelope/bremsstrahlung/pdebr" << Z <<
".p08";
389 ost << path <<
"/penelope/bremsstrahlung/pdebr0" << Z <<
".p08";
390 std::ifstream
file(ost.str().c_str());
393 G4String excep =
"G4PenelopeBremsstrahlungFS - data file " +
394 G4String(ost.str()) +
" not found!";
395 G4Exception(
"G4PenelopeBremsstrahlungFS::ReadDataFile()",
407 ed <<
"Corrupted data file for Z=" << Z <<
G4endl;
408 G4Exception(
"G4PenelopeBremsstrahlungFS::ReadDataFile()",
415 for (
size_t ie=0;ie<
nBinsE;ie++)
422 for (
size_t ix=0;ix<
nBinsX;ix++)
425 (*theMatrix)[ie*(nBinsX+1)+ix] = myDouble*
millibarn;
428 (*theMatrix)[ie*(nBinsX+1)+nBinsX] = myDouble*
millibarn;
449 size_t size = nBinsX;
453 if (momOrder<-1 || size<2 || theXGrid[0]<0)
455 G4Exception(
"G4PenelopeBremsstrahlungFS::GetMomentumIntegral()",
459 for (
size_t i=1;i<size;i++)
461 if (theXGrid[i]<0 || theXGrid[i]<theXGrid[i-1])
464 ed <<
"Invalid call for bin " << i <<
G4endl;
465 G4Exception(
"G4PenelopeBremsstrahlungFS::GetMomentumIntegral()",
472 if (xup < theXGrid[0])
474 bool loopAgain =
true;
477 for (
size_t i=0;i<size-1;i++)
493 if (std::fabs(dx)>1
e-14*std::fabs(dy))
498 ds = a*std::log(xtc/x1)+b*(xtc-
x1);
499 else if (momOrder == 0)
500 ds = a*(xtc-
x1) + 0.5*b*(xtc*xtc-x1*x1);
502 ds = a*(std::pow(xtc,momOrder+1)-std::pow(x1,momOrder+1))/((
G4double) (momOrder + 1))
503 + b*(std::pow(xtc,momOrder+2)-std::pow(x1,momOrder+2))/((
G4double) (momOrder + 2));
506 ds = 0.5*(y1+
y2)*(xtc-x1)*std::pow(xtc,momOrder);
520 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
524 G4Exception(
"G4PenelopeBremsstrahlungFS::GetScaledXSTable()",
525 "em2013",
FatalException,
"Unable to retrieve the cross section table");
537 G4cout <<
"Entering in G4PenelopeBremsstrahlungFS::InitializeEnergySampling() for " <<
541 std::pair<const G4Material*,G4double> theKey = std::make_pair(material,cut);
549 for (
size_t i=0;i<
nBinsE;i++)
555 G4Exception(
"G4PenelopeBremsstrahlungFS::InitializeEnergySampling()",
556 "em2013",
FatalException,
"Unable to retrieve the cross section table");
559 for (
size_t ie=0;ie<
nBinsE;ie++)
566 for (
size_t ix=1;ix<
nBinsX;ix++)
590 for (
size_t ix=0;ix<
nBinsX;ix++)
593 tempData[ix] =
G4Exp((*vv)[ie+1]);
603 thePBcut->insert(std::make_pair(theKey,thePBvec));
612 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
616 ed <<
"Unable to retrieve the SamplingTable: " <<
619 G4Exception(
"G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
629 G4bool firstOrLastBin =
false;
634 firstOrLastBin =
true;
639 firstOrLastBin =
true;
672 G4cout <<
"Creating new instance of G4PhysicsFreeVector() on the worker" <<
G4endl;
680 for (
size_t iloop=0;iloop<
nBinsX;iloop++)
682 G4double val = (*theVec1)[iloop]+(((*theVec2)[iloop]-(*theVec1)[iloop]))*
689 for (
size_t iloop=0;iloop<
nBinsX;iloop++)
698 pbcut = (*(
thePBcut->find(theKey)->second))[eBin] +
699 ((*(
thePBcut->find(theKey)->second))[eBin+1]-(*(
thePBcut->find(theKey)->second))[eBin])*
706 G4int nIterations = 0;
714 if (pt < (*theTempVec)[0])
716 else if (pt > (*theTempVec)[
nBinsX-1])
721 if (delta < pt*1
e-10)
725 ed <<
"Found that (pt > (*theTempVec)[nBinsX-1]) with pt = " << pt <<
726 " , (*theTempVec)[nBinsX-1] = " << (*theTempVec)[
nBinsX-1] <<
" and delta = " <<
728 ed <<
"Possible symptom of problem with numerical precision" <<
G4endl;
729 G4Exception(
"G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
735 ed <<
"Crash at (pt > (*theTempVec)[nBinsX-1]) with pt = " << pt <<
736 " , (*theTempVec)[nBinsX-1]=" << (*theTempVec)[
nBinsX-1] <<
" and nBinsX = " <<
738 ed <<
"Material: " << mat->
GetName() <<
", energy = " << energy/
keV <<
" keV" <<
740 G4Exception(
"G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
751 if (pt > (*theTempVec)[
k])
789 ed <<
"Warning in G4PenelopeBremsstrahlungFS::SampleX()" <<
G4endl;
790 ed <<
"Conflicting end-point values: w1=" << w1 <<
"; w2 = " << w2 <<
G4endl;
791 ed <<
"wbcut = " << wbcut <<
" energy= " << energy/
keV <<
" keV" <<
G4endl;
792 ed <<
"cut = " << cut/
keV <<
" keV" <<
G4endl;
793 G4Exception(
"G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
"em2015",
809 if (nIterations > 100)
811 }
while(eGamma < cut);