51 theReducedXSTable(0),theEffectiveZSq(0),theSamplingTable(0),
55 {1.0e-12,0.025e0,0.05e0,0.075e0,0.1e0,0.15e0,0.2e0,0.25e0,
56 0.3e0,0.35e0,0.4e0,0.45e0,0.5e0,0.55e0,0.6e0,0.65e0,0.7e0,
57 0.75e0,0.8e0,0.85e0,0.9e0,0.925e0,0.95e0,0.97e0,0.99e0,
58 0.995e0,0.999e0,0.9995e0,0.9999e0,0.99995e0,0.99999e0,1.0e0};
60 for (
size_t ix=0;ix<nBinsX;ix++)
61 theXGrid[ix] = tempvector[ix];
63 for (
size_t i=0;i<nBinsE;i++)
66 theElementData =
new std::map<G4int,G4DataVector*>;
80 std::map<G4int,G4DataVector*>::iterator i;
83 for (i=theElementData->begin(); i != theElementData->end(); i++)
85 delete theElementData;
96 std::map< std::pair<const G4Material*,G4double> ,
G4PhysicsTable*>::iterator j;
98 if (theReducedXSTable)
100 for (j=theReducedXSTable->begin(); j != theReducedXSTable->end(); j++)
106 delete theReducedXSTable;
107 theReducedXSTable = 0;
110 if (theSamplingTable)
112 for (j=theSamplingTable->begin(); j != theSamplingTable->end(); j++)
118 delete theSamplingTable;
119 theSamplingTable = 0;
125 for (kk=thePBcut->begin(); kk != thePBcut->end(); kk++)
134 delete theEffectiveZSq;
145 if (!theEffectiveZSq)
148 ed <<
"The container for the <Z^2> values is not initialized" <<
G4endl;
149 G4Exception(
"G4PenelopeBremsstrahlungFS::GetEffectiveZSquared()",
154 if (theEffectiveZSq->count(material))
155 return theEffectiveZSq->find(material)->second;
159 ed <<
"The value of <Z^2> is not properly set for material " <<
162 G4Exception(
"G4PenelopeBremsstrahlungFS::GetEffectiveZSquared()",
170void G4PenelopeBremsstrahlungFS::BuildScaledXSTable(
const G4Material* material,
183 std::vector<G4double> *StechiometricFactors =
new std::vector<G4double>;
187 for (
G4int i=0;i<nElements;i++)
189 G4double fraction = fractionVector[i];
190 G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole);
191 StechiometricFactors->push_back(fraction/atomicWeigth);
194 G4double MaxStechiometricFactor = 0.;
195 for (
G4int i=0;i<nElements;i++)
197 if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
198 MaxStechiometricFactor = (*StechiometricFactors)[i];
201 for (
G4int i=0;i<nElements;i++)
202 (*StechiometricFactors)[i] /= MaxStechiometricFactor;
206 for (
G4int i=0;i<nElements;i++)
208 G4double Z = (*elementVector)[i]->GetZ();
209 sumz2 += (*StechiometricFactors)[i]*Z*Z;
210 sums += (*StechiometricFactors)[i];
214 theEffectiveZSq->insert(std::make_pair(material,ZBR2));
223 for (
G4int iel=0;iel<nElements;iel++)
225 G4double Z = (*elementVector)[iel]->GetZ();
227 G4double wgt = (*StechiometricFactors)[iel]*Z*Z/ZBR2;
231 if (!theElementData->count(iZ))
234 if (!theElementData->count(iZ))
237 ed <<
"Error in G4PenelopeBremsstrahlungFS::BuildScaledXSTable" <<
G4endl;
238 ed <<
"Unable to retrieve data for element " << iZ <<
G4endl;
239 G4Exception(
"G4PenelopeBremsstrahlungFS::BuildScaledXSTable()",
244 G4DataVector* atomData = theElementData->find(iZ)->second;
246 for (
size_t ie=0;ie<nBinsE;ie++)
248 (*tempData)[ie] += wgt*(*atomData)[ie*(nBinsX+1)+nBinsX];
249 for (
size_t ix=0;ix<nBinsX;ix++)
250 (*tempMatrix)[ie*nBinsX+ix] += wgt*(*atomData)[ie*(nBinsX+1)+ix];
258 for (
size_t ie=0;ie<nBinsE;ie++)
262 for (
size_t ix=0;ix<nBinsX;ix++)
263 tempData2[ix] = (*tempMatrix)[ie*nBinsX+ix];
266 G4double fact = millibarn*(theEGrid[ie]+electron_mass_c2)*(1./fine_structure_const)/
267 (classic_electr_radius*classic_electr_radius*(theEGrid[ie]+2.0*electron_mass_c2));
268 G4double fnorm = (*tempData)[ie]/(rsum*fact);
269 G4double TST = 100.*std::fabs(fnorm-1.0);
273 ed <<
"G4PenelopeBremsstrahlungFS. Corrupted data files?" <<
G4endl;
274 G4cout <<
"TST= " << TST <<
"; fnorm = " << fnorm <<
G4endl;
277 G4cout << ie <<
" " << theEGrid[ie]/keV <<
" " << (*tempData)[ie]/barn <<
G4endl;
278 G4Exception(
"G4PenelopeBremsstrahlungFS::BuildScaledXSTable()",
281 for (
size_t ix=0;ix<nBinsX;ix++)
282 (*tempMatrix)[ie*nBinsX+ix] *= fnorm;
296 for (
size_t i=0;i<nBinsX;i++)
299 for (
size_t ix=0;ix<nBinsX;ix++)
303 for (
size_t ie=0;ie<nBinsE;ie++)
305 G4double logene = std::log(theEGrid[ie]);
306 G4double aValue = (*tempMatrix)[ie*nBinsX+ix];
307 if (aValue < 1e-20*millibarn)
308 aValue = 1e-20*millibarn;
309 theVec->
PutValue(ie+1,logene,std::log(aValue));
315 G4double val1eV = (*theVec)[1]+derivative*(log1eV-theVec->
Energy(1));
319 std::pair<const G4Material*,G4double> theKey = std::make_pair(material,cut);
320 theReducedXSTable->insert(std::make_pair(theKey,thePhysicsTable));
322 delete StechiometricFactors;
331void G4PenelopeBremsstrahlungFS::ReadDataFile(
G4int Z)
334 char* path = getenv(
"G4LEDATA");
337 G4String excep =
"G4PenelopeBremsstrahlungFS - G4LEDATA environment variable not set!";
338 G4Exception(
"G4PenelopeBremsstrahlungFS::ReadDataFile()",
345 std::ostringstream ost;
347 ost << path <<
"/penelope/bremsstrahlung/pdebr" << Z <<
".p08";
349 ost << path <<
"/penelope/bremsstrahlung/pdebr0" << Z <<
".p08";
350 std::ifstream file(ost.str().c_str());
353 G4String excep =
"G4PenelopeBremsstrahlungFS - data file " +
354 G4String(ost.str()) +
" not found!";
355 G4Exception(
"G4PenelopeBremsstrahlungFS::ReadDataFile()",
367 ed <<
"Corrupted data file for Z=" << Z <<
G4endl;
368 G4Exception(
"G4PenelopeBremsstrahlungFS::ReadDataFile()",
375 for (
size_t ie=0;ie<nBinsE;ie++)
380 theEGrid[ie] = myDouble*eV;
382 for (
size_t ix=0;ix<nBinsX;ix++)
385 (*theMatrix)[ie*(nBinsX+1)+ix] = myDouble*millibarn;
388 (*theMatrix)[ie*(nBinsX+1)+nBinsX] = myDouble*millibarn;
392 theElementData->insert(std::make_pair(Z,theMatrix));
409 size_t size = nBinsX;
413 if (momOrder<-1 || size<2 || theXGrid[0]<0)
415 G4Exception(
"G4PenelopeBremsstrahlungFS::GetMomentumIntegral()",
419 for (
size_t i=1;i<size;i++)
421 if (theXGrid[i]<0 || theXGrid[i]<theXGrid[i-1])
424 ed <<
"Invalid call for bin " << i <<
G4endl;
425 G4Exception(
"G4PenelopeBremsstrahlungFS::GetMomentumIntegral()",
432 if (xup < theXGrid[0])
434 bool loopAgain =
true;
435 G4double xt = std::min(xup,theXGrid[size-1]);
437 for (
size_t i=0;i<size-1;i++)
439 G4double x1 = std::max(theXGrid[i],eps);
441 G4double x2 = std::max(theXGrid[i+1],eps);
453 if (std::fabs(dx)>1e-14*std::fabs(dy))
458 ds = a*std::log(xtc/x1)+b*(xtc-x1);
459 else if (momOrder == 0)
460 ds = a*(xtc-x1) + 0.5*b*(xtc*xtc-x1*x1);
462 ds = a*(std::pow(xtc,momOrder+1)-std::pow(x1,momOrder+1))/((
G4double) (momOrder + 1))
463 + b*(std::pow(xtc,momOrder+2)-std::pow(x1,momOrder+2))/((
G4double) (momOrder + 2));
466 ds = 0.5*(y1+y2)*(xtc-x1)*std::pow(xtc,momOrder);
480 if (!theReducedXSTable)
481 theReducedXSTable =
new std::map< std::pair<const G4Material*,G4double> ,
483 if (!theEffectiveZSq)
484 theEffectiveZSq =
new std::map<const G4Material*,G4double>;
487 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
488 if (!(theReducedXSTable->count(theKey)))
489 BuildScaledXSTable(mat,cut);
491 if (!(theReducedXSTable->count(theKey)))
493 G4Exception(
"G4PenelopeBremsstrahlungFS::GetScaledXSTable()",
494 "em2013",
FatalException,
"Unable to retrieve the cross section table");
497 return theReducedXSTable->find(theKey)->second;
502void G4PenelopeBremsstrahlungFS::InitializeEnergySampling(
const G4Material* material,
505 std::pair<const G4Material*,G4double> theKey = std::make_pair(material,cut);
513 for (
size_t i=0;i<nBinsE;i++)
521 for (
size_t ie=0;ie<nBinsE;ie++)
527 theVec->
PutValue(0,theXGrid[0],value);
528 for (
size_t ix=1;ix<nBinsX;ix++)
535 G4double x1=std::max(theXGrid[ix-1],1.0e-35);
539 G4double x2=std::max(theXGrid[ix],1.0e-35);
543 G4double dS = A*std::log(x2/x1)+B*(x2-x1);
545 theVec->
PutValue(ix,theXGrid[ix],value);
552 for (
size_t ix=0;ix<nBinsX;ix++)
555 tempData[ix] = std::exp((*vv)[ie+1]);
560 thePBvec->
PutValue(ie,theEGrid[ie],pbval);
564 theSamplingTable->insert(std::make_pair(theKey,thePhysicsTable));
565 thePBcut->insert(std::make_pair(theKey,thePBvec));
574 if (!theSamplingTable)
576 new std::map< std::pair<const G4Material*,G4double> ,
G4PhysicsTable*>;
581 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
583 if (!(theSamplingTable->count(theKey)))
585 InitializeEnergySampling(mat,cut);
586 if (!(theSamplingTable->count(theKey)) || !(thePBcut->count(theKey)))
589 ed <<
"Unable to create the SamplingTable: " <<
590 theSamplingTable->count(theKey) <<
" " <<
591 thePBcut->count(theKey) <<
G4endl;
592 G4Exception(
"G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
597 G4PhysicsTable* theTableInte = theSamplingTable->find(theKey)->second;
598 G4PhysicsTable* theTableRed = theReducedXSTable->find(theKey)->second;
602 G4bool firstOrLastBin =
false;
604 if (energy < theEGrid[0])
607 firstOrLastBin =
true;
609 else if (energy > theEGrid[nBinsE-1])
612 firstOrLastBin =
true;
621 if (energy > theEGrid[k])
640 for (
size_t iloop=0;iloop<nBinsX;iloop++)
642 G4double val = (*theVec1)[iloop]+(((*theVec2)[iloop]-(*theVec1)[iloop]))*
643 (energy-theEGrid[eBin])/(theEGrid[eBin+1]-theEGrid[eBin]);
644 theTempVec->
PutValue(iloop,theXGrid[iloop],val);
649 for (
size_t iloop=0;iloop<nBinsX;iloop++)
650 theTempVec->
PutValue(iloop,theXGrid[iloop],(*theVec1)[iloop]);
654 G4double pbcut = (*(thePBcut->find(theKey)->second))[eBin];
658 pbcut = (*(thePBcut->find(theKey)->second))[eBin] +
659 ((*(thePBcut->find(theKey)->second))[eBin+1]-(*(thePBcut->find(theKey)->second))[eBin])*
660 (energy-theEGrid[eBin])/(theEGrid[eBin+1]-theEGrid[eBin]);
663 G4double pCumulative = (*theTempVec)[nBinsX-1];
672 if (pt < (*theTempVec)[0])
674 else if (pt > (*theTempVec)[nBinsX-1])
678 G4double delta = pt-(*theTempVec)[nBinsX-1];
679 if (delta < pt*1e-10)
683 ed <<
"Found that (pt > (*theTempVec)[nBinsX-1]) with pt = " << pt <<
684 " , (*theTempVec)[nBinsX-1] = " << (*theTempVec)[nBinsX-1] <<
" and delta = " <<
685 (pt-(*theTempVec)[nBinsX-1]) <<
G4endl;
686 ed <<
"Possible symptom of problem with numerical precision" <<
G4endl;
687 G4Exception(
"G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
693 ed <<
"Crash at (pt > (*theTempVec)[nBinsX-1]) with pt = " << pt <<
694 " , (*theTempVec)[nBinsX-1]=" << (*theTempVec)[nBinsX-1] <<
" and nBinsX = " <<
696 ed <<
"Material: " << mat->
GetName() <<
", energy = " << energy/keV <<
" keV" <<
698 G4Exception(
"G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
709 if (pt > (*theTempVec)[k])
724 G4double pdf1 = std::exp((*v1)[eBin+1]);
725 G4double pdf2 = std::exp((*v2)[eBin+1]);
732 G4double wbcut = (cut < energy) ? cut/energy : 1.0;
734 wbcut = (cut < theEGrid[eBin]) ? cut/theEGrid[eBin] : 1.0;
740 G4cout <<
"Warning in G4PenelopeBremsstrahlungFS::SampleX()" <<
G4endl;
741 G4cout <<
"Conflicting end-point values: w1=" << w1 <<
"; w2 = " << w2 <<
G4endl;
742 G4cout <<
"wbcut = " << wbcut <<
" energy= " << energy/keV <<
" keV" <<
G4endl;
747 G4double pmax = std::max(A+B*w1,A+B*w2);
757 }
while(eGamma < cut);
std::vector< G4Element * > G4ElementVector
G4DLLIMPORT std::ostream G4cout
const G4ElementVector * GetElementVector() const
const G4double * GetFractionVector() const
size_t GetNumberOfElements() const
const G4String & GetName() const
G4double GetMomentumIntegral(G4double *y, G4double up, G4int momOrder)
G4double SampleGammaEnergy(G4double energy, const G4Material *, G4double cut)
G4double GetEffectiveZSquared(const G4Material *mat)
~G4PenelopeBremsstrahlungFS()
G4PhysicsTable * GetScaledXSTable(const G4Material *, G4double cut)
G4PenelopeBremsstrahlungFS()
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
void push_back(G4PhysicsVector *)
G4double Energy(size_t index) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
std::ostringstream G4ExceptionDescription