G4PenelopeBremsstrahlungModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // Author: Luciano Pandola
00029 //
00030 // History:
00031 // --------
00032 // 23 Nov 2010   L Pandola    First complete implementation, Penelope v2008
00033 // 24 May 2011   L. Pandola   Renamed (make default Penelope)
00034 // 13 Mar 2012   L. Pandola   Updated the interface for the angular generator
00035 // 18 Jul 2012   L. Pandola   Migrate to the new interface of the angular generator, which 
00036 //                            now provides the G4ThreeVector and takes care of rotation
00037 //
00038 
00039 #include "G4PenelopeBremsstrahlungModel.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042 #include "G4PenelopeBremsstrahlungFS.hh"
00043 #include "G4PenelopeBremsstrahlungAngular.hh"
00044 #include "G4ParticleDefinition.hh"
00045 #include "G4MaterialCutsCouple.hh"
00046 #include "G4ProductionCutsTable.hh"
00047 #include "G4DynamicParticle.hh"
00048 #include "G4Gamma.hh"
00049 #include "G4Electron.hh"
00050 #include "G4Positron.hh"
00051 #include "G4PenelopeOscillatorManager.hh"
00052 #include "G4PenelopeCrossSection.hh"
00053 #include "G4PhysicsFreeVector.hh"
00054 #include "G4PhysicsLogVector.hh" 
00055 #include "G4PhysicsTable.hh"
00056 
00057 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00058  
00059 G4PenelopeBremsstrahlungModel::G4PenelopeBremsstrahlungModel(const G4ParticleDefinition*,
00060                                                              const G4String& nam)
00061   :G4VEmModel(nam),fParticleChange(0),isInitialised(false),energyGrid(0),  
00062    XSTableElectron(0),XSTablePositron(0)
00063   
00064 {
00065   fIntrinsicLowEnergyLimit = 100.0*eV;
00066   fIntrinsicHighEnergyLimit = 100.0*GeV;
00067   nBins = 200;
00068 
00069   SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
00070   //
00071   oscManager = G4PenelopeOscillatorManager::GetOscillatorManager();
00072   //
00073   verboseLevel= 0;
00074    
00075   // Verbosity scale:
00076   // 0 = nothing
00077   // 1 = warning for energy non-conservation
00078   // 2 = details of energy budget
00079   // 3 = calculation of cross sections, file openings, sampling of atoms
00080   // 4 = entering in methods
00081 
00082   fPenelopeFSHelper = new G4PenelopeBremsstrahlungFS();
00083   fPenelopeAngular = new G4PenelopeBremsstrahlungAngular();
00084 
00085   // Atomic deexcitation model activated by default
00086   SetDeexcitationFlag(true);
00087 
00088 }
00089 
00090 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00091  
00092 G4PenelopeBremsstrahlungModel::~G4PenelopeBremsstrahlungModel()
00093 {
00094   ClearTables();
00095   if (fPenelopeFSHelper)
00096     delete fPenelopeFSHelper;
00097   if (fPenelopeAngular)
00098     delete fPenelopeAngular;
00099 }
00100 
00101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00102 
00103 void G4PenelopeBremsstrahlungModel::Initialise(const G4ParticleDefinition*,
00104                                              const G4DataVector&)
00105 {
00106   if (verboseLevel > 3)
00107     G4cout << "Calling G4PenelopeBremsstrahlungModel::Initialise()" << G4endl;
00108  
00109   //Clear and re-build the tables 
00110   ClearTables();
00111   
00112   //forces the cleaning of tables, in this specific case
00113   if (fPenelopeAngular)
00114     fPenelopeAngular->Initialize();
00115     
00116 
00117   //Set the number of bins for the tables. 20 points per decade
00118   nBins = (size_t) (20*std::log10(HighEnergyLimit()/LowEnergyLimit()));
00119   nBins = std::max(nBins,(size_t)100);
00120   energyGrid = new G4PhysicsLogVector(LowEnergyLimit(),
00121                                       HighEnergyLimit(), 
00122                                       nBins-1); //one hidden bin is added
00123  
00124 
00125   XSTableElectron = new 
00126     std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
00127   XSTablePositron = new 
00128     std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
00129 
00130   if (verboseLevel > 2) {
00131     G4cout << "Penelope Bremsstrahlung model v2008 is initialized " << G4endl
00132            << "Energy range: "
00133            << LowEnergyLimit() / keV << " keV - "
00134            << HighEnergyLimit() / GeV << " GeV." 
00135            << G4endl;
00136   }
00137  
00138   if(isInitialised) return;
00139   fParticleChange = GetParticleChangeForLoss();
00140   isInitialised = true;
00141 }
00142 
00143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00144 
00145 G4double G4PenelopeBremsstrahlungModel::CrossSectionPerVolume(const G4Material* material,
00146                                            const G4ParticleDefinition* theParticle,
00147                                            G4double energy,
00148                                            G4double cutEnergy,
00149                                            G4double)
00150 {
00151   //
00152   if (verboseLevel > 3)
00153     G4cout << "Calling CrossSectionPerVolume() of G4PenelopeBremsstrahlungModel" << G4endl;
00154  
00155   SetupForMaterial(theParticle, material, energy);
00156 
00157   G4double crossPerMolecule = 0.;
00158 
00159   G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(theParticle,material,
00160                                                                 cutEnergy);
00161 
00162   if (theXS)
00163     crossPerMolecule = theXS->GetHardCrossSection(energy);
00164 
00165   G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
00166   G4double atPerMol =  oscManager->GetAtomsPerMolecule(material);
00167  
00168   if (verboseLevel > 3)
00169     G4cout << "Material " << material->GetName() << " has " << atPerMol <<
00170       "atoms per molecule" << G4endl;
00171 
00172   G4double moleculeDensity = 0.; 
00173   if (atPerMol)
00174     moleculeDensity = atomDensity/atPerMol;
00175  
00176   G4double crossPerVolume = crossPerMolecule*moleculeDensity;
00177 
00178   if (verboseLevel > 2)
00179   {
00180     G4cout << "G4PenelopeBremsstrahlungModel " << G4endl;
00181     G4cout << "Mean free path for gamma emission > " << cutEnergy/keV << " keV at " <<
00182       energy/keV << " keV = " << (1./crossPerVolume)/mm << " mm" << G4endl;
00183   }
00184   
00185   return crossPerVolume;
00186 }
00187 
00188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00189                                                                                                      
00190 //This is a dummy method. Never inkoved by the tracking, it just issues
00191 //a warning if one tries to get Cross Sections per Atom via the
00192 //G4EmCalculator.
00193 G4double G4PenelopeBremsstrahlungModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
00194                                                                      G4double,
00195                                                                      G4double,
00196                                                                      G4double,
00197                                                                      G4double,
00198                                                                      G4double)
00199 {
00200   G4cout << "*** G4PenelopeBremsstrahlungModel -- WARNING ***" << G4endl;
00201   G4cout << "Penelope Bremsstrahlung model v2008 does not calculate cross section _per atom_ " << G4endl;
00202   G4cout << "so the result is always zero. For physics values, please invoke " << G4endl;
00203   G4cout << "GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" << G4endl;
00204   return 0;
00205 }
00206 
00207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00208 
00209 G4double G4PenelopeBremsstrahlungModel::ComputeDEDXPerVolume(const G4Material* material,
00210                                                              const G4ParticleDefinition* theParticle,
00211                                                              G4double kineticEnergy,
00212                                                              G4double cutEnergy)
00213 {
00214   if (verboseLevel > 3)
00215     G4cout << "Calling ComputeDEDX() of G4PenelopeBremsstrahlungModel" << G4endl;
00216  
00217   G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(theParticle,material,
00218                                                                 cutEnergy);
00219   G4double sPowerPerMolecule = 0.0;
00220   if (theXS)
00221     sPowerPerMolecule = theXS->GetSoftStoppingPower(kineticEnergy);
00222 
00223   
00224   G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
00225   G4double atPerMol =  oscManager->GetAtomsPerMolecule(material);
00226                                                                                         
00227   G4double moleculeDensity = 0.; 
00228   if (atPerMol)
00229     moleculeDensity = atomDensity/atPerMol;
00230  
00231   G4double sPowerPerVolume = sPowerPerMolecule*moleculeDensity;
00232 
00233   if (verboseLevel > 2)
00234     {
00235       G4cout << "G4PenelopeBremsstrahlungModel " << G4endl;
00236       G4cout << "Stopping power < " << cutEnergy/keV << " keV at " <<
00237         kineticEnergy/keV << " keV = " << 
00238         sPowerPerVolume/(keV/mm) << " keV/mm" << G4endl;
00239     }
00240   return sPowerPerVolume;
00241 }
00242 
00243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00244 
00245 void G4PenelopeBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>*fvect,
00246                                                       const G4MaterialCutsCouple* couple,
00247                                                       const G4DynamicParticle* aDynamicParticle,
00248                                                       G4double cutG,
00249                                                       G4double)
00250 {
00251   if (verboseLevel > 3)
00252     G4cout << "Calling SampleSecondaries() of G4PenelopeBremsstrahlungModel" << G4endl;
00253 
00254   G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
00255   const G4Material* material = couple->GetMaterial();
00256 
00257   if (kineticEnergy <= fIntrinsicLowEnergyLimit)
00258    {
00259      fParticleChange->SetProposedKineticEnergy(0.);
00260      fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
00261      return ;
00262    }
00263 
00264   G4ParticleMomentum particleDirection0 = aDynamicParticle->GetMomentumDirection();
00265   //This is the momentum
00266   G4ThreeVector initialMomentum =  aDynamicParticle->GetMomentum();
00267 
00268   //Not enough energy to produce a secondary! Return with nothing happened
00269   if (kineticEnergy < cutG)
00270     return;
00271 
00272   if (verboseLevel > 3)
00273     G4cout << "Going to sample gamma energy for: " <<material->GetName() << " " << 
00274       "energy = " << kineticEnergy/keV << ", cut = " << cutG/keV << G4endl;
00275 
00276    //Sample gamma's energy according to the spectrum
00277   G4double gammaEnergy = 
00278     fPenelopeFSHelper->SampleGammaEnergy(kineticEnergy,material,cutG);
00279 
00280   if (verboseLevel > 3)
00281     G4cout << "Sampled gamma energy: " << gammaEnergy/keV << " keV" << G4endl;
00282   
00283   //Now sample the direction for the Gamma. Notice that the rotation is already done
00284   //Z is unused here, I plug 0. The information is in the material pointer
00285    G4ThreeVector gammaDirection1 = 
00286      fPenelopeAngular->SampleDirection(aDynamicParticle,gammaEnergy,0,material);
00287   
00288   if (verboseLevel > 3)
00289     G4cout << "Sampled cosTheta for e-: " << gammaDirection1.cosTheta() << G4endl;
00290 
00291   G4double residualPrimaryEnergy = kineticEnergy-gammaEnergy;
00292   if (residualPrimaryEnergy < 0)
00293     {
00294       //Ok we have a problem, all energy goes with the gamma
00295       gammaEnergy += residualPrimaryEnergy;
00296       residualPrimaryEnergy = 0.0;
00297     }
00298   
00299   //Produce final state according to momentum conservation
00300   G4ThreeVector particleDirection1 = initialMomentum - gammaEnergy*gammaDirection1;
00301   particleDirection1 = particleDirection1.unit(); //normalize
00302 
00303   //Update the primary particle
00304   if (residualPrimaryEnergy > 0.)
00305     {
00306       fParticleChange->ProposeMomentumDirection(particleDirection1);
00307       fParticleChange->SetProposedKineticEnergy(residualPrimaryEnergy);
00308     }
00309   else
00310     {
00311       fParticleChange->SetProposedKineticEnergy(0.);
00312     }
00313   
00314   //Now produce the photon
00315   G4DynamicParticle* theGamma = new G4DynamicParticle(G4Gamma::Gamma(),
00316                                                       gammaDirection1,
00317                                                       gammaEnergy);
00318   fvect->push_back(theGamma);
00319   
00320   if (verboseLevel > 1)
00321     {
00322       G4cout << "-----------------------------------------------------------" << G4endl;
00323       G4cout << "Energy balance from G4PenelopeBremsstrahlung" << G4endl;
00324       G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
00325       G4cout << "-----------------------------------------------------------" << G4endl;
00326       G4cout << "Outgoing primary energy: " << residualPrimaryEnergy/keV << " keV" << G4endl;
00327       G4cout << "Bremsstrahlung photon " << gammaEnergy/keV << " keV" << G4endl;
00328       G4cout << "Total final state: " << (residualPrimaryEnergy+gammaEnergy)/keV 
00329              << " keV" << G4endl;
00330       G4cout << "-----------------------------------------------------------" << G4endl;
00331     }
00332 
00333   if (verboseLevel > 0)
00334     {
00335       G4double energyDiff = std::fabs(residualPrimaryEnergy+gammaEnergy-kineticEnergy);
00336       if (energyDiff > 0.05*keV)
00337         G4cout << "Warning from G4PenelopeBremsstrahlung: problem with energy conservation: " 
00338                <<
00339           (residualPrimaryEnergy+gammaEnergy)/keV <<
00340           " keV (final) vs. " <<
00341           kineticEnergy/keV << " keV (initial)" << G4endl;
00342     }  
00343   return;
00344 }
00345 
00346 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00347 
00348 void G4PenelopeBremsstrahlungModel::ClearTables()
00349 {  
00350   std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>::iterator i;
00351   if (XSTableElectron)
00352     {
00353       for (i=XSTableElectron->begin(); i != XSTableElectron->end(); i++)
00354         {
00355           G4PenelopeCrossSection* tab = i->second;
00356           delete tab;
00357         }
00358       delete XSTableElectron;
00359       XSTableElectron = 0;
00360     }
00361   if (XSTablePositron)
00362     {
00363       for (i=XSTablePositron->begin(); i != XSTablePositron->end(); i++)
00364         {
00365           G4PenelopeCrossSection* tab = i->second;
00366           delete tab;
00367         }
00368       delete XSTablePositron;
00369       XSTablePositron = 0;
00370     } 
00371 
00372   if (energyGrid)
00373     delete energyGrid;
00374 
00375   if (fPenelopeFSHelper)
00376     fPenelopeFSHelper->ClearTables(); 
00377 
00378   if (verboseLevel > 2)
00379     G4cout << "G4PenelopeBremsstrahlungModel: cleared tables" << G4endl;
00380 
00381  return;
00382 }
00383 
00384 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00385 
00386 G4double G4PenelopeBremsstrahlungModel::MinEnergyCut(const G4ParticleDefinition*,
00387                                                        const G4MaterialCutsCouple*)
00388 {
00389   return fIntrinsicLowEnergyLimit;
00390 }
00391 
00392 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00393 
00394 void G4PenelopeBremsstrahlungModel::BuildXSTable(const G4Material* mat,G4double cut)
00395 {
00396 //
00397   //This method fills the G4PenelopeCrossSection containers for electrons or positrons
00398   //and for the given material/cut couple. 
00399   //Equivalent of subroutines EBRaT and PINaT of Penelope
00400   //
00401   if (verboseLevel > 2)
00402     {
00403       G4cout << "G4PenelopeBremsstrahlungModel: going to build cross section table " << G4endl;
00404       G4cout << "for e+/e- in " << mat->GetName() << G4endl;
00405     }
00406 
00407   //Tables have been already created (checked by GetCrossSectionTableForCouple)
00408   if (energyGrid->GetVectorLength() != nBins) 
00409     {
00410       G4ExceptionDescription ed;
00411       ed << "Energy Grid looks not initialized" << G4endl;
00412       ed << nBins << " " << energyGrid->GetVectorLength() << G4endl;
00413       G4Exception("G4PenelopeBremsstrahlungModel::BuildXSTable()",
00414                   "em2016",FatalException,ed);
00415     }
00416 
00417   G4PenelopeCrossSection* XSEntry = new G4PenelopeCrossSection(nBins);
00418   G4PenelopeCrossSection* XSEntry2 = new G4PenelopeCrossSection(nBins);
00419 
00420   G4PhysicsTable* table = fPenelopeFSHelper->GetScaledXSTable(mat,cut);
00421 
00422 
00423   //loop on the energy grid  
00424   for (size_t bin=0;bin<nBins;bin++)
00425     {
00426        G4double energy = energyGrid->GetLowEdgeEnergy(bin);
00427        G4double XH0=0, XH1=0, XH2=0;
00428        G4double XS0=0, XS1=0, XS2=0;
00429     
00430        //Global xs factor
00431        G4double fact = fPenelopeFSHelper->GetEffectiveZSquared(mat)*
00432          ((energy+electron_mass_c2)*(energy+electron_mass_c2)/
00433           (energy*(energy+2.0*electron_mass_c2))); 
00434        
00435        G4double restrictedCut = cut/energy;
00436    
00437        //Now I need the dSigma/dX profile - interpolated on energy - for 
00438        //the 32-point x grid. Interpolation is log-log
00439        size_t nBinsX = fPenelopeFSHelper->GetNBinsX();       
00440        G4double* tempData = new G4double[nBinsX];     
00441        G4double logene = std::log(energy);
00442        for (size_t ix=0;ix<nBinsX;ix++) 
00443          {
00444            //find dSigma/dx for the given E. X belongs to the 32-point grid.
00445            G4double val = (*table)[ix]->Value(logene);     
00446            tempData[ix] = std::exp(val); //back to the real value!
00447          } 
00448        
00449        G4double XH0A = 0.;
00450        if (restrictedCut <= 1) //calculate only if we are above threshold!
00451          XH0A = fPenelopeFSHelper->GetMomentumIntegral(tempData,1.0,-1) -
00452            fPenelopeFSHelper->GetMomentumIntegral(tempData,restrictedCut,-1);
00453        G4double XS1A = fPenelopeFSHelper->GetMomentumIntegral(tempData,
00454                                                               restrictedCut,0);       
00455        G4double XS2A = fPenelopeFSHelper->GetMomentumIntegral(tempData,
00456                                                               restrictedCut,1);
00457        G4double XH1A=0, XH2A=0;
00458        if (restrictedCut <=1)
00459          {
00460            XH1A = fPenelopeFSHelper->GetMomentumIntegral(tempData,1.0,0) -
00461              XS1A;
00462            XH2A = fPenelopeFSHelper->GetMomentumIntegral(tempData,1.0,1) -
00463              XS2A;
00464          }
00465        delete[] tempData;
00466 
00467        XH0 = XH0A*fact;
00468        XS1 = XS1A*fact*energy;
00469        XH1 = XH1A*fact*energy;
00470        XS2 = XS2A*fact*energy*energy;
00471        XH2 = XH2A*fact*energy*energy;
00472 
00473        XSEntry->AddCrossSectionPoint(bin,energy,XH0,XH1,XH2,XS0,XS1,XS2);
00474 
00475        //take care of positrons
00476        G4double posCorrection = GetPositronXSCorrection(mat,energy);
00477        XSEntry2->AddCrossSectionPoint(bin,energy,XH0*posCorrection,
00478                                       XH1*posCorrection,
00479                                       XH2*posCorrection,
00480                                       XS0,
00481                                       XS1*posCorrection,
00482                                       XS2*posCorrection);
00483     }
00484  
00485   //Insert in the appropriate table
00486   std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);  
00487   XSTableElectron->insert(std::make_pair(theKey,XSEntry));
00488   XSTablePositron->insert(std::make_pair(theKey,XSEntry2));
00489   
00490   return;
00491 }
00492 
00493 
00494 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00495 
00496 G4PenelopeCrossSection* 
00497 G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple(const G4ParticleDefinition* part,
00498                                                              const G4Material* mat,
00499                                                              G4double cut)
00500 {
00501   if (part != G4Electron::Electron() && part != G4Positron::Positron())
00502     {
00503       G4ExceptionDescription ed;
00504       ed << "Invalid particle: " << part->GetParticleName() << G4endl;
00505       G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
00506                   "em0001",FatalException,ed);
00507       return NULL;
00508     }
00509 
00510   if (part == G4Electron::Electron())
00511     {
00512       if (!XSTableElectron)
00513         {         
00514           G4String excep = "The Cross Section Table for e- was not initialized correctly!";
00515           G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
00516                       "em2013",FatalException,excep);          
00517           return NULL;
00518         }
00519       std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
00520       if (XSTableElectron->count(theKey)) //table already built
00521         return XSTableElectron->find(theKey)->second;
00522       else
00523         {
00524           BuildXSTable(mat,cut);
00525           if (XSTableElectron->count(theKey)) //now it should be ok!
00526             return XSTableElectron->find(theKey)->second;
00527           else
00528             {
00529               G4ExceptionDescription ed;
00530               ed << "Unable to build e- table for " << mat->GetName() << G4endl;
00531               G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
00532                           "em2009",FatalException,ed);
00533             }
00534         }
00535     }
00536   if (part == G4Positron::Positron())
00537     {
00538       if (!XSTablePositron)
00539         {
00540           G4String excep = "The Cross Section Table for e+ was not initialized correctly!";
00541           G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
00542                       "em2013",FatalException,excep); 
00543           return NULL;
00544         }
00545       std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
00546       if (XSTablePositron->count(theKey)) //table already built
00547         return XSTablePositron->find(theKey)->second;
00548       else
00549         {
00550           BuildXSTable(mat,cut);
00551           if (XSTablePositron->count(theKey)) //now it should be ok!
00552             return XSTablePositron->find(theKey)->second;
00553           else
00554             {
00555               G4ExceptionDescription ed;
00556               ed << "Unable to build e+ table for " << mat->GetName() << G4endl;
00557               G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
00558                           "em2009",FatalException,ed);
00559             }
00560         }
00561     }
00562   return NULL;
00563 }
00564 
00565 
00566 
00567 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00568 
00569 G4double G4PenelopeBremsstrahlungModel::GetPositronXSCorrection(const G4Material* mat,
00570                                                                   G4double energy)
00571 {
00572   //The electron-to-positron correction factor is set equal to the ratio of the 
00573   //radiative stopping powers for positrons and electrons, which has been calculated 
00574   //by Kim et al. (1986) (cf. Berger and Seltzer, 1982). Here, it is used an 
00575   //analytical approximation which reproduces the tabulated values with 0.5% 
00576   //accuracy
00577   G4double t=std::log(1.0+1e6*energy/
00578                       (electron_mass_c2*fPenelopeFSHelper->GetEffectiveZSquared(mat)));
00579   G4double corr = 1.0-std::exp(-t*(1.2359e-1-t*(6.1274e-2-t*
00580                                            (3.1516e-2-t*(7.7446e-3-t*(1.0595e-3-t*
00581                                                                       (7.0568e-5-t*
00582                                                                        1.8080e-6)))))));
00583   return corr;
00584 }

Generated on Mon May 27 17:49:17 2013 for Geant4 by  doxygen 1.4.7