G4SeltzerBergerModel Class Reference

#include <G4SeltzerBergerModel.hh>

Inheritance diagram for G4SeltzerBergerModel:

G4eBremsstrahlungRelModel G4VEmModel

Public Member Functions

 G4SeltzerBergerModel (const G4ParticleDefinition *p=0, const G4String &nam="eBremSB")
virtual ~G4SeltzerBergerModel ()
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double cutEnergy, G4double maxEnergy)
void SetBicubicInterpolationFlag (G4bool)

Protected Member Functions

virtual G4double ComputeDXSectionPerAtom (G4double gammaEnergy)

Detailed Description

Definition at line 61 of file G4SeltzerBergerModel.hh.

Constructor & Destructor Documentation

G4SeltzerBergerModel::G4SeltzerBergerModel ( const G4ParticleDefinition *  p = 0,
const G4String &  nam = "eBremSB" 
)

Definition at line 78 of file G4SeltzerBergerModel.cc.

References G4VEmModel::SetLowEnergyLimit(), and G4VEmModel::SetLPMFlag().

00080   : G4eBremsstrahlungRelModel(p,nam),useBicubicInterpolation(false)
00081 {
00082   SetLowEnergyLimit(0.0);
00083   SetLPMFlag(false);
00084   nwarn = 0;
00085 }

G4SeltzerBergerModel::~G4SeltzerBergerModel (  )  [virtual]

Definition at line 89 of file G4SeltzerBergerModel.cc.

00090 {
00091   for(size_t i=0; i<101; ++i) { 
00092     if(dataSB[i]) {
00093       delete dataSB[i]; 
00094       dataSB[i] = 0;
00095     } 
00096   }
00097 }

Member Function Documentation

G4double G4SeltzerBergerModel::ComputeDXSectionPerAtom ( G4double  gammaEnergy  )  [protected, virtual]

Reimplemented from G4eBremsstrahlungRelModel.

Definition at line 175 of file G4SeltzerBergerModel.cc.

References G4eBremsstrahlungRelModel::bremFactor, G4eBremsstrahlungRelModel::currentZ, G4eBremsstrahlungRelModel::isElectron, G4eBremsstrahlungRelModel::kinEnergy, G4eBremsstrahlungRelModel::particleMass, G4eBremsstrahlungRelModel::totalEnergy, and G4Physics2DVector::Value().

00176 {
00177 
00178   if(gammaEnergy < 0.0 || kinEnergy <= 0.0) { return 0.0; }
00179   G4double x = gammaEnergy/kinEnergy;
00180   G4double y = log(kinEnergy/MeV);
00181   G4int Z = G4int(currentZ);
00182   //G4cout << "G4SeltzerBergerModel::ComputeDXSectionPerAtom Z= " << Z
00183   //     << " x= " << x << " y= " << y << " " << dataSB[Z] << G4endl;
00184   if(!dataSB[Z]) { ReadData(Z); }
00185   G4double invb2 = totalEnergy*totalEnergy/(kinEnergy*(kinEnergy + 2*particleMass));
00186   G4double cross = dataSB[Z]->Value(x,y)*invb2*millibarn/bremFactor;
00187   
00188   if(!isElectron) {
00189     G4double invbeta1 = sqrt(invb2);
00190     G4double e2 = kinEnergy - gammaEnergy;
00191     if(e2 > 0.0) {
00192       G4double invbeta2 = (e2 + particleMass)/sqrt(e2*(e2 + 2*particleMass));
00193       G4double xxx = twopi*fine_structure_const*currentZ*(invbeta1 - invbeta2);
00194       if(xxx < expnumlim) { cross = 0.0; }
00195       else { cross *= exp(xxx); }
00196     } else {
00197       cross = 0.0;
00198     }
00199   }
00200   
00201   return cross;
00202 }

void G4SeltzerBergerModel::Initialise ( const G4ParticleDefinition *  ,
const G4DataVector &   
) [virtual]

Reimplemented from G4eBremsstrahlungRelModel.

Definition at line 101 of file G4SeltzerBergerModel.cc.

References G4Element::GetElementTable(), G4Element::GetNumberOfElements(), and G4eBremsstrahlungRelModel::Initialise().

00103 {
00104   // check environment variable
00105   // Build the complete string identifying the file with the data set
00106   char* path = getenv("G4LEDATA");
00107 
00108   // Access to elements
00109   const G4ElementTable* theElmTable = G4Element::GetElementTable();
00110   size_t numOfElm = G4Element::GetNumberOfElements();
00111   if(numOfElm > 0) {
00112     for(size_t i=0; i<numOfElm; ++i) {
00113       G4int Z = G4int(((*theElmTable)[i])->GetZ());
00114       if(Z < 1)        { Z = 1; }
00115       else if(Z > 100) { Z = 100; }
00116       //G4cout << "Z= " << Z << G4endl;
00117       // Initialisation
00118       if(!dataSB[Z]) { ReadData(Z, path); }
00119     }
00120   }
00121 
00122   G4eBremsstrahlungRelModel::Initialise(p, cuts);
00123 }

void G4SeltzerBergerModel::SampleSecondaries ( std::vector< G4DynamicParticle * > *  ,
const G4MaterialCutsCouple *  ,
const G4DynamicParticle *  ,
G4double  cutEnergy,
G4double  maxEnergy 
) [virtual]

Reimplemented from G4eBremsstrahlungRelModel.

Definition at line 207 of file G4SeltzerBergerModel.cc.

References G4eBremsstrahlungRelModel::currentZ, G4eBremsstrahlungRelModel::densityCorr, G4eBremsstrahlungRelModel::densityFactor, G4eBremsstrahlungRelModel::fParticleChange, fStopAndKill, G4cout, G4endl, G4UniformRand, G4VEmModel::GetAngularDistribution(), G4DynamicParticle::GetKineticEnergy(), G4MaterialCutsCouple::GetMaterial(), G4DynamicParticle::GetMomentumDirection(), G4ParticleDefinition::GetParticleName(), G4eBremsstrahlungRelModel::isElectron, G4eBremsstrahlungRelModel::particle, G4eBremsstrahlungRelModel::particleMass, G4VParticleChange::ProposeTrackStatus(), G4VEmAngularDistribution::SampleDirection(), G4VEmModel::SecondaryThreshold(), G4VEmModel::SelectRandomAtom(), G4eBremsstrahlungRelModel::SetCurrentElement(), G4ParticleChangeForLoss::SetProposedKineticEnergy(), G4ParticleChangeForLoss::SetProposedMomentumDirection(), G4eBremsstrahlungRelModel::SetupForMaterial(), G4eBremsstrahlungRelModel::theGamma, G4eBremsstrahlungRelModel::totalEnergy, G4VEmModel::Value(), and G4Physics2DVector::Value().

00212 {
00213   G4double kineticEnergy = dp->GetKineticEnergy();
00214   G4double cut  = std::min(cutEnergy, kineticEnergy);
00215   G4double emax = std::min(maxEnergy, kineticEnergy);
00216   if(cut >= emax) { return; }
00217 
00218   SetupForMaterial(particle, couple->GetMaterial(), kineticEnergy);
00219 
00220   const G4Element* elm = 
00221     SelectRandomAtom(couple,particle,kineticEnergy,cut,emax);
00222   SetCurrentElement(elm->GetZ());
00223   G4int Z = G4int(currentZ);
00224 
00225   totalEnergy = kineticEnergy + particleMass;
00226   densityCorr = densityFactor*totalEnergy*totalEnergy;
00227   G4double totMomentum = sqrt(kineticEnergy*(totalEnergy + electron_mass_c2));
00228   /*
00229   G4cout << "G4SeltzerBergerModel::SampleSecondaries E(MeV)= " 
00230          << kineticEnergy/MeV
00231          << " Z= " << Z << " cut(MeV)= " << cut/MeV 
00232          << " emax(MeV)= " << emax/MeV << " corr= " << densityCorr << G4endl;
00233   */
00234   G4double xmin = log(cut*cut + densityCorr);
00235   G4double xmax = log(emax*emax  + densityCorr);
00236   G4double y = log(kineticEnergy/MeV);
00237 
00238   G4double gammaEnergy, v; 
00239 
00240   // majoranta
00241   G4double x0 = cut/kineticEnergy;
00242   G4double vmax = dataSB[Z]->Value(x0, y)*1.02;
00243   //  G4double invbeta1 = 0;
00244 
00245   const G4double epeaklimit= 300*MeV; 
00246   const G4double elowlimit = 10*keV; 
00247 
00248   // majoranta corrected for e-
00249   if(isElectron && x0 < 0.97 && 
00250      ((kineticEnergy > epeaklimit) || (kineticEnergy < elowlimit))) {
00251     G4double ylim = std::min(ylimit[Z],1.1*dataSB[Z]->Value(0.97, y)); 
00252     if(ylim > vmax) { vmax = ylim; }
00253   }
00254   if(x0 < 0.05) { vmax *= 1.2; }
00255 
00256   //G4cout<<"y= "<<y<<" xmin= "<<xmin<<" xmax= "<<xmax<<" vmax= "<<vmax<<G4endl;
00257   //  G4int ncount = 0;
00258   do {
00259     //++ncount;
00260     G4double x = exp(xmin + G4UniformRand()*(xmax - xmin)) - densityCorr;
00261     if(x < 0.0) { x = 0.0; }
00262     gammaEnergy = sqrt(x);
00263     G4double x1 = gammaEnergy/kineticEnergy;
00264     v = dataSB[Z]->Value(x1, y);
00265 
00266     // correction for positrons        
00267     if(!isElectron) {
00268       G4double e1 = kineticEnergy - cut;
00269       G4double invbeta1 = (e1 + particleMass)/sqrt(e1*(e1 + 2*particleMass));
00270       G4double e2 = kineticEnergy - gammaEnergy;
00271       G4double invbeta2 = (e2 + particleMass)/sqrt(e2*(e2 + 2*particleMass));
00272       G4double xxx = twopi*fine_structure_const*currentZ*(invbeta1 - invbeta2); 
00273 
00274       if(xxx < expnumlim) { v = 0.0; }
00275       else { v *= exp(xxx); }
00276     }
00277    
00278     if (v > 1.05*vmax && nwarn < 20) {
00279       ++nwarn;
00280       G4cout << "### G4SeltzerBergerModel Warning: Majoranta exceeded! "
00281              << v << " > " << vmax << " by " << v/vmax
00282              << " Egamma(MeV)= " << gammaEnergy
00283              << " Ee(MeV)= " << kineticEnergy
00284              << " Z= " << Z << "  " << particle->GetParticleName()
00285         //<< " ncount= " << ncount
00286              << G4endl;
00287      
00288       if ( 20 == nwarn ) {
00289         G4cout << "### G4SeltzerBergerModel Warnings will not be printed anymore"
00290                << G4endl;
00291       }
00292     }
00293   } while (v < vmax*G4UniformRand());
00294 
00295   //
00296   // angles of the emitted gamma. ( Z - axis along the parent particle)
00297   // use general interface
00298   //
00299 
00300   G4ThreeVector gammaDirection = 
00301     GetAngularDistribution()->SampleDirection(dp, totalEnergy-gammaEnergy,
00302                                               Z, couple->GetMaterial());
00303 
00304   // create G4DynamicParticle object for the Gamma
00305   G4DynamicParticle* gamma = 
00306     new G4DynamicParticle(theGamma,gammaDirection,gammaEnergy);
00307   vdp->push_back(gamma);
00308   
00309   G4ThreeVector direction = (totMomentum*dp->GetMomentumDirection()
00310                              - gammaEnergy*gammaDirection).unit();
00311 
00312   /*
00313   G4cout << "### G4SBModel: v= "
00314          << " Eg(MeV)= " << gammaEnergy
00315          << " Ee(MeV)= " << kineticEnergy
00316          << " DirE " << direction << " DirG " << gammaDirection
00317          << G4endl;
00318   */
00319   // energy of primary
00320   G4double finalE = kineticEnergy - gammaEnergy;
00321 
00322   // stop tracking and create new secondary instead of primary
00323   if(gammaEnergy > SecondaryThreshold()) {
00324     fParticleChange->ProposeTrackStatus(fStopAndKill);
00325     fParticleChange->SetProposedKineticEnergy(0.0);
00326     G4DynamicParticle* el = 
00327       new G4DynamicParticle(const_cast<G4ParticleDefinition*>(particle),
00328                             direction, finalE);
00329     vdp->push_back(el);
00330 
00331     // continue tracking
00332   } else {
00333     fParticleChange->SetProposedMomentumDirection(direction);
00334     fParticleChange->SetProposedKineticEnergy(finalE);
00335   }
00336 }

void G4SeltzerBergerModel::SetBicubicInterpolationFlag ( G4bool   )  [inline]

Definition at line 100 of file G4SeltzerBergerModel.hh.

00101 {
00102   useBicubicInterpolation = val;
00103 }

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