G4PenelopeBremsstrahlungFS Class Reference

#include <G4PenelopeBremsstrahlungFS.hh>


Public Member Functions

 G4PenelopeBremsstrahlungFS ()
 ~G4PenelopeBremsstrahlungFS ()
G4double GetEffectiveZSquared (const G4Material *mat)
void ClearTables ()
size_t GetNBinsX ()
G4double GetMomentumIntegral (G4double *y, G4double up, G4int momOrder)
G4PhysicsTable * GetScaledXSTable (const G4Material *, G4double cut)
G4double SampleGammaEnergy (G4double energy, const G4Material *, G4double cut)

Detailed Description

Definition at line 55 of file G4PenelopeBremsstrahlungFS.hh.

Constructor & Destructor Documentation

G4PenelopeBremsstrahlungFS::G4PenelopeBremsstrahlungFS (  ) 

Definition at line 50 of file G4PenelopeBremsstrahlungFS.cc.

00050                                                        : 
00051   theReducedXSTable(0),theEffectiveZSq(0),theSamplingTable(0),
00052   thePBcut(0)
00053 {
00054   G4double tempvector[nBinsX] = 
00055     {1.0e-12,0.025e0,0.05e0,0.075e0,0.1e0,0.15e0,0.2e0,0.25e0,
00056     0.3e0,0.35e0,0.4e0,0.45e0,0.5e0,0.55e0,0.6e0,0.65e0,0.7e0,
00057     0.75e0,0.8e0,0.85e0,0.9e0,0.925e0,0.95e0,0.97e0,0.99e0,
00058     0.995e0,0.999e0,0.9995e0,0.9999e0,0.99995e0,0.99999e0,1.0e0};
00059 
00060   for (size_t ix=0;ix<nBinsX;ix++)
00061     theXGrid[ix] = tempvector[ix];
00062 
00063   for (size_t i=0;i<nBinsE;i++)
00064     theEGrid[i] = 0.;
00065 
00066   theElementData = new std::map<G4int,G4DataVector*>;
00067   theTempVec = new G4PhysicsFreeVector(nBinsX);
00068 }

G4PenelopeBremsstrahlungFS::~G4PenelopeBremsstrahlungFS (  ) 

Definition at line 72 of file G4PenelopeBremsstrahlungFS.cc.

References ClearTables().

00073 {
00074   ClearTables();
00075 
00076   if (theTempVec)
00077     delete theTempVec;
00078 
00079   //Clear manually theElementData
00080   std::map<G4int,G4DataVector*>::iterator i;
00081   if (theElementData)
00082     {
00083       for (i=theElementData->begin(); i != theElementData->end(); i++)        
00084         delete i->second;        
00085       delete theElementData;
00086       theElementData = 0;
00087     }
00088 
00089 }

Member Function Documentation

void G4PenelopeBremsstrahlungFS::ClearTables (  ) 

Definition at line 94 of file G4PenelopeBremsstrahlungFS.cc.

Referenced by ~G4PenelopeBremsstrahlungFS().

00095 {  
00096   std::map< std::pair<const G4Material*,G4double> ,G4PhysicsTable*>::iterator j;
00097 
00098   if (theReducedXSTable)
00099     {
00100       for (j=theReducedXSTable->begin(); j != theReducedXSTable->end(); j++)
00101         {
00102           G4PhysicsTable* tab = j->second;
00103           tab->clearAndDestroy();
00104           delete tab;
00105         }
00106       delete theReducedXSTable;
00107       theReducedXSTable = 0;
00108     }
00109 
00110   if (theSamplingTable)
00111     {
00112       for (j=theSamplingTable->begin(); j != theSamplingTable->end(); j++)
00113         {
00114           G4PhysicsTable* tab = j->second;
00115           tab->clearAndDestroy();
00116           delete tab;
00117         }
00118       delete theSamplingTable;
00119       theSamplingTable = 0;
00120     }
00121 
00122   std::map< std::pair<const G4Material*,G4double> ,G4PhysicsFreeVector*>::iterator kk;
00123   if (thePBcut)
00124     {
00125       for (kk=thePBcut->begin(); kk != thePBcut->end(); kk++)   
00126         delete kk->second;
00127       delete thePBcut;
00128       thePBcut = 0;      
00129     }
00130 
00131 
00132   if (theEffectiveZSq)
00133     {
00134       delete theEffectiveZSq;
00135       theEffectiveZSq = 0;
00136     }
00137   
00138  return;
00139 }

G4double G4PenelopeBremsstrahlungFS::GetEffectiveZSquared ( const G4Material *  mat  ) 

Definition at line 143 of file G4PenelopeBremsstrahlungFS.cc.

References FatalException, G4endl, and G4Exception().

00144 {
00145   if (!theEffectiveZSq)
00146     {
00147       G4ExceptionDescription ed;
00148       ed << "The container for the <Z^2> values is not initialized" << G4endl;
00149       G4Exception("G4PenelopeBremsstrahlungFS::GetEffectiveZSquared()",
00150                   "em2007",FatalException,ed);
00151       return 0;
00152     }
00153   //found in the table: return it 
00154   if (theEffectiveZSq->count(material))
00155     return theEffectiveZSq->find(material)->second;    
00156   else
00157     {
00158       G4ExceptionDescription ed;
00159       ed << "The value of  <Z^2> is not properly set for material " << 
00160         material->GetName() << G4endl;
00161       //requires running of BuildScaledXSTable()
00162       G4Exception("G4PenelopeBremsstrahlungFS::GetEffectiveZSquared()",
00163                   "em2008",FatalException,ed);
00164     }
00165   return 0;
00166 }

G4double G4PenelopeBremsstrahlungFS::GetMomentumIntegral ( G4double *  y,
G4double  up,
G4int  momOrder 
)

Definition at line 400 of file G4PenelopeBremsstrahlungFS.cc.

References FatalException, G4endl, and G4Exception().

00403 {
00404   //Corresponds to the function RLMOM of Penelope
00405   //This method performs the calculation of the integral of (x^momOrder)*y over the interval
00406   //from x[0] to xup, obtained by linear interpolation on a table of y. 
00407   //The independent variable is assumed to take positive values only.
00408   //
00409   size_t size = nBinsX;
00410   const G4double eps = 1e-35;
00411 
00412   //Check that the call is valid
00413   if (momOrder<-1 || size<2 || theXGrid[0]<0)
00414     {
00415       G4Exception("G4PenelopeBremsstrahlungFS::GetMomentumIntegral()",
00416                   "em2011",FatalException,"Invalid call");
00417     }
00418 
00419   for (size_t i=1;i<size;i++)
00420     {
00421       if (theXGrid[i]<0 || theXGrid[i]<theXGrid[i-1])
00422         {
00423           G4ExceptionDescription ed;
00424           ed << "Invalid call for bin " << i << G4endl;
00425           G4Exception("G4PenelopeBremsstrahlungFS::GetMomentumIntegral()",
00426                   "em2012",FatalException,ed);
00427         }
00428     }
00429 
00430   //Compute the integral
00431   G4double result = 0;
00432   if (xup < theXGrid[0])
00433     return result;
00434   bool loopAgain = true;
00435   G4double xt = std::min(xup,theXGrid[size-1]);
00436   G4double xtc = 0;
00437   for (size_t i=0;i<size-1;i++)
00438     {
00439       G4double x1 = std::max(theXGrid[i],eps);
00440       G4double y1 = y[i];
00441       G4double x2 = std::max(theXGrid[i+1],eps);
00442       G4double y2 = y[i+1];
00443       if (xt < x2)
00444         {
00445           xtc = xt;
00446           loopAgain = false;
00447         }
00448       else
00449         xtc = x2;
00450       G4double dx = x2-x1;
00451       G4double dy = y2-y1;
00452       G4double ds = 0;
00453       if (std::fabs(dx)>1e-14*std::fabs(dy))
00454         {
00455           G4double b=dy/dx;
00456           G4double a=y1-b*x1;    
00457           if (momOrder == -1)       
00458             ds = a*std::log(xtc/x1)+b*(xtc-x1);
00459           else if (momOrder == 0) //speed it up, not using pow()
00460             ds = a*(xtc-x1) + 0.5*b*(xtc*xtc-x1*x1);
00461           else
00462             ds = a*(std::pow(xtc,momOrder+1)-std::pow(x1,momOrder+1))/((G4double) (momOrder + 1))
00463               + b*(std::pow(xtc,momOrder+2)-std::pow(x1,momOrder+2))/((G4double) (momOrder + 2));           
00464         }
00465       else
00466         ds = 0.5*(y1+y2)*(xtc-x1)*std::pow(xtc,momOrder);
00467       result += ds;
00468       if (!loopAgain) 
00469         return result;
00470     }
00471   return result;
00472 }

size_t G4PenelopeBremsstrahlungFS::GetNBinsX (  )  [inline]

Definition at line 64 of file G4PenelopeBremsstrahlungFS.hh.

00064 {return nBinsX;};

G4PhysicsTable * G4PenelopeBremsstrahlungFS::GetScaledXSTable ( const G4Material *  ,
G4double  cut 
)

Definition at line 476 of file G4PenelopeBremsstrahlungFS.cc.

References FatalException, and G4Exception().

00478 {
00479   //check if the container exists (if not, create it)
00480   if (!theReducedXSTable)
00481     theReducedXSTable = new std::map< std::pair<const G4Material*,G4double> ,
00482     G4PhysicsTable*>;
00483   if (!theEffectiveZSq)
00484     theEffectiveZSq = new std::map<const G4Material*,G4double>;
00485 
00486   //check if it already contains the entry
00487   std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
00488   if (!(theReducedXSTable->count(theKey))) //not found
00489     BuildScaledXSTable(mat,cut);
00490   
00491   if (!(theReducedXSTable->count(theKey)))
00492     {
00493       G4Exception("G4PenelopeBremsstrahlungFS::GetScaledXSTable()",
00494                   "em2013",FatalException,"Unable to retrieve the cross section table");
00495     }
00496 
00497   return theReducedXSTable->find(theKey)->second;
00498 }

G4double G4PenelopeBremsstrahlungFS::SampleGammaEnergy ( G4double  energy,
const G4Material *  ,
G4double  cut 
)

Definition at line 571 of file G4PenelopeBremsstrahlungFS.cc.

References FatalException, G4cout, G4endl, G4Exception(), G4UniformRand, G4Material::GetName(), JustWarning, and G4PhysicsFreeVector::PutValue().

Referenced by G4PenelopeBremsstrahlungModel::SampleSecondaries().

00573 {
00574   if (!theSamplingTable)    
00575     theSamplingTable = 
00576       new std::map< std::pair<const G4Material*,G4double> , G4PhysicsTable*>;
00577   if (!thePBcut)
00578     thePBcut = 
00579       new std::map< std::pair<const G4Material*,G4double> , G4PhysicsFreeVector* >;
00580 
00581   std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
00582   
00583   if (!(theSamplingTable->count(theKey)))
00584     {
00585       InitializeEnergySampling(mat,cut);
00586       if (!(theSamplingTable->count(theKey)) || !(thePBcut->count(theKey)))
00587         {
00588           G4ExceptionDescription ed;
00589           ed << "Unable to create the SamplingTable: " << 
00590             theSamplingTable->count(theKey) << " " << 
00591             thePBcut->count(theKey) << G4endl;
00592           G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
00593                       "em2014",FatalException,ed);        
00594         }
00595     }
00596 
00597   G4PhysicsTable* theTableInte = theSamplingTable->find(theKey)->second;
00598   G4PhysicsTable* theTableRed = theReducedXSTable->find(theKey)->second;
00599 
00600   //Find the energy bin using bi-partition
00601   size_t eBin = 0;
00602   G4bool firstOrLastBin = false;
00603 
00604   if (energy < theEGrid[0]) //below first bin
00605     {
00606       eBin = 0;
00607       firstOrLastBin = true;
00608     }
00609   else if (energy > theEGrid[nBinsE-1]) //after last bin
00610     {
00611       eBin = nBinsE-1;
00612       firstOrLastBin = true;
00613     }
00614   else
00615     {
00616       size_t i=0;
00617       size_t j=nBinsE-1;
00618       while ((j-i)>1)
00619         {
00620           size_t k = (i+j)/2;
00621           if (energy > theEGrid[k])
00622             i = k;
00623           else
00624             j = k;       
00625         }
00626       eBin = i;
00627     }
00628 
00629   //Get the appropriate physics vector
00630   G4PhysicsFreeVector* theVec1 = (G4PhysicsFreeVector*) (*theTableInte)[eBin];
00631 
00632   //use a "temporary" vector which contains the linear interpolation of the x spectra 
00633   //in energy
00634 
00635   //theTempVect is allocated only once (member variable), but it is overwritten at 
00636   //every call of this method (because the interpolation factors change!)
00637   if (!firstOrLastBin)
00638     {  
00639       G4PhysicsFreeVector* theVec2 = (G4PhysicsFreeVector*) (*theTableInte)[eBin+1];
00640       for (size_t iloop=0;iloop<nBinsX;iloop++)
00641         {
00642           G4double val = (*theVec1)[iloop]+(((*theVec2)[iloop]-(*theVec1)[iloop]))*
00643             (energy-theEGrid[eBin])/(theEGrid[eBin+1]-theEGrid[eBin]);
00644           theTempVec->PutValue(iloop,theXGrid[iloop],val);
00645         }      
00646     }
00647   else //first or last bin, no interpolation
00648     {
00649       for (size_t iloop=0;iloop<nBinsX;iloop++) 
00650         theTempVec->PutValue(iloop,theXGrid[iloop],(*theVec1)[iloop]);  
00651     }
00652 
00653   //Start the game
00654   G4double pbcut = (*(thePBcut->find(theKey)->second))[eBin];
00655  
00656   if (!firstOrLastBin) //linear interpolation on pbcut as well
00657     {
00658       pbcut = (*(thePBcut->find(theKey)->second))[eBin] + 
00659         ((*(thePBcut->find(theKey)->second))[eBin+1]-(*(thePBcut->find(theKey)->second))[eBin])*
00660         (energy-theEGrid[eBin])/(theEGrid[eBin+1]-theEGrid[eBin]);
00661     }
00662 
00663   G4double pCumulative = (*theTempVec)[nBinsX-1]; //last value  
00664  
00665   G4double eGamma = 0;
00666   do
00667     {    
00668       G4double pt = pbcut + G4UniformRand()*(pCumulative - pbcut);
00669 
00670       //find where it is
00671       size_t ibin = 0;
00672       if (pt < (*theTempVec)[0])
00673         ibin = 0;
00674       else if (pt > (*theTempVec)[nBinsX-1])
00675         {
00676           //We observed problems due to numerical rounding here (STT). 
00677           //delta here is a tiny positive number
00678           G4double delta = pt-(*theTempVec)[nBinsX-1];
00679           if (delta < pt*1e-10) // very small! Numerical rounding only
00680             {
00681               ibin = nBinsX-2;
00682               G4ExceptionDescription ed;
00683               ed << "Found that (pt > (*theTempVec)[nBinsX-1]) with pt = " << pt << 
00684                 " , (*theTempVec)[nBinsX-1] = " << (*theTempVec)[nBinsX-1] << " and delta = " << 
00685                 (pt-(*theTempVec)[nBinsX-1]) << G4endl;
00686               ed << "Possible symptom of problem with numerical precision" << G4endl;
00687               G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
00688                           "em2015",JustWarning,ed);
00689             }
00690           else //real problem
00691             {
00692               G4ExceptionDescription ed;
00693               ed << "Crash at (pt > (*theTempVec)[nBinsX-1]) with pt = " << pt << 
00694                 " , (*theTempVec)[nBinsX-1]=" << (*theTempVec)[nBinsX-1] << " and nBinsX = " << 
00695                 nBinsX << G4endl;
00696               ed << "Material: " << mat->GetName() << ", energy = " << energy/keV << " keV" << 
00697                 G4endl;
00698               G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
00699                       "em2015",FatalException,ed);        
00700             }
00701         }
00702       else
00703         {
00704           size_t i=0;
00705           size_t j=nBinsX-1;
00706           while ((j-i)>1)
00707             {
00708               size_t k = (i+j)/2;
00709               if (pt > (*theTempVec)[k])
00710                 i = k;
00711               else
00712                 j = k;
00713             }
00714           ibin = i;
00715         }
00716     
00717       G4double w1 = theXGrid[ibin];
00718       G4double w2 = theXGrid[ibin+1];            
00719 
00720       G4PhysicsFreeVector* v1 = (G4PhysicsFreeVector*) (*theTableRed)[ibin];
00721       G4PhysicsFreeVector* v2 = (G4PhysicsFreeVector*) (*theTableRed)[ibin+1];     
00722       //Remember: the table theReducedXSTable has a fake first point in energy
00723       //so, it contains one more bin than nBinsE.
00724       G4double pdf1 = std::exp((*v1)[eBin+1]);
00725       G4double pdf2 = std::exp((*v2)[eBin+1]);    
00726       G4double deltaW = w2-w1;
00727       G4double dpdfb = pdf2-pdf1;
00728       G4double B = dpdfb/deltaW;
00729       G4double A = pdf1-B*w1;
00730       //I already made an interpolation in energy, so I can use the actual value for the 
00731       //calculation of the wbcut, instead of the grid values (except for the last bin)
00732       G4double wbcut  = (cut < energy) ? cut/energy : 1.0;
00733       if (firstOrLastBin) //this is an particular case: no interpolation available
00734         wbcut  = (cut < theEGrid[eBin]) ? cut/theEGrid[eBin] : 1.0;
00735      
00736       if (w1 < wbcut)
00737         w1 = wbcut;
00738       if (w2 < w1)
00739         {
00740           G4cout << "Warning in G4PenelopeBremsstrahlungFS::SampleX()" << G4endl;
00741           G4cout << "Conflicting end-point values: w1=" << w1 << "; w2 = " << w2 << G4endl;
00742           G4cout << "wbcut = " << wbcut << " energy= " << energy/keV << " keV" << G4endl;
00743           G4cout << "cut = " << cut/keV << " keV" << G4endl;
00744           return w1*energy;
00745         }
00746   
00747       G4double pmax = std::max(A+B*w1,A+B*w2);
00748       G4bool loopAgain = false;    
00749       do
00750         {
00751           loopAgain = false;
00752           eGamma = w1* std::pow((w2/w1),G4UniformRand());
00753           if  (G4UniformRand()*pmax > (A+B*eGamma))
00754             loopAgain = true;   
00755         }while(loopAgain);     
00756       eGamma *= energy;
00757     }while(eGamma < cut); //repeat if sampled sub-cut!  
00758 
00759   return eGamma;
00760 }

The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:52 2013 for Geant4 by  doxygen 1.4.7