G4Nucleus.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 //
00027 //
00028  // original by H.P. Wellisch
00029  // modified by J.L. Chuma, TRIUMF, 19-Nov-1996
00030  // last modified: 27-Mar-1997
00031  // J.P.Wellisch: 23-Apr-97: minor simplifications
00032  // modified by J.L.Chuma 24-Jul-97  to set the total momentum in Cinema and
00033  //                                  EvaporationEffects
00034  // modified by J.L.Chuma 21-Oct-97  put std::abs() around the totalE^2-mass^2
00035  //                                  in calculation of total momentum in
00036  //                                  Cinema and EvaporationEffects
00037  // Chr. Volcker, 10-Nov-1997: new methods and class variables.
00038  // HPW added utilities for low energy neutron transport. (12.04.1998)
00039  // M.G. Pia, 2 Oct 1998: modified GetFermiMomentum to avoid memory leaks
00040  // G.Folger, spring 2010:  add integer A/Z interface
00041  
00042 #include "G4Nucleus.hh"
00043 #include "G4NucleiProperties.hh"
00044 #include "G4PhysicalConstants.hh"
00045 #include "G4SystemOfUnits.hh"
00046 #include "Randomize.hh"
00047 #include "G4HadronicException.hh"
00048  
00049 G4Nucleus::G4Nucleus()
00050   : theA(0), theZ(0), aEff(0.0), zEff(0)
00051 {
00052   pnBlackTrackEnergy = 0.0;
00053   dtaBlackTrackEnergy = 0.0;
00054   pnBlackTrackEnergyfromAnnihilation = 0.0;
00055   dtaBlackTrackEnergyfromAnnihilation = 0.0;
00056   excitationEnergy = 0.0;
00057   momentum = G4ThreeVector(0.,0.,0.);
00058   fermiMomentum = 1.52*hbarc/fermi;
00059   theTemp = 293.16*kelvin;
00060   fIsotope = 0;
00061 }
00062 
00063 G4Nucleus::G4Nucleus( const G4double A, const G4double Z )
00064 {
00065   SetParameters( A, Z );
00066   pnBlackTrackEnergy = 0.0;
00067   dtaBlackTrackEnergy = 0.0;
00068   pnBlackTrackEnergyfromAnnihilation = 0.0;
00069   dtaBlackTrackEnergyfromAnnihilation = 0.0;
00070   excitationEnergy = 0.0;
00071   momentum = G4ThreeVector(0.,0.,0.);
00072   fermiMomentum = 1.52*hbarc/fermi;
00073   theTemp = 293.16*kelvin;
00074   fIsotope = 0;
00075 }
00076 
00077 G4Nucleus::G4Nucleus( const G4int A, const G4int Z )
00078 {
00079   SetParameters( A, Z );
00080   pnBlackTrackEnergy = 0.0;
00081   dtaBlackTrackEnergy = 0.0;
00082   pnBlackTrackEnergyfromAnnihilation = 0.0;
00083   dtaBlackTrackEnergyfromAnnihilation = 0.0;
00084   excitationEnergy = 0.0;
00085   momentum = G4ThreeVector(0.,0.,0.);
00086   fermiMomentum = 1.52*hbarc/fermi;
00087   theTemp = 293.16*kelvin;
00088   fIsotope = 0;
00089 }
00090 
00091 G4Nucleus::G4Nucleus( const G4Material *aMaterial )
00092 {
00093   ChooseParameters( aMaterial );
00094   pnBlackTrackEnergy = 0.0;
00095   dtaBlackTrackEnergy = 0.0;
00096   pnBlackTrackEnergyfromAnnihilation = 0.0;
00097   dtaBlackTrackEnergyfromAnnihilation = 0.0;
00098   excitationEnergy = 0.0;
00099   momentum = G4ThreeVector(0.,0.,0.);
00100   fermiMomentum = 1.52*hbarc/fermi;
00101   theTemp = aMaterial->GetTemperature();
00102   fIsotope = 0;
00103 }
00104 
00105 G4Nucleus::~G4Nucleus() {}
00106 
00107 G4ReactionProduct G4Nucleus::
00108 GetBiasedThermalNucleus(G4double aMass, G4ThreeVector aVelocity, G4double temp) const
00109 {
00110   G4double velMag = aVelocity.mag();
00111   G4ReactionProduct result;
00112   G4double value = 0;
00113   G4double random = 1;
00114   G4double norm = 3.*std::sqrt(k_Boltzmann*temp*aMass*G4Neutron::Neutron()->GetPDGMass());
00115   norm /= G4Neutron::Neutron()->GetPDGMass();
00116   norm *= 5.;
00117   norm += velMag;
00118   norm /= velMag;
00119   while(value/norm<random)
00120   {
00121      result = GetThermalNucleus(aMass, temp);
00122      G4ThreeVector targetVelocity = 1./result.GetMass()*result.GetMomentum();
00123      value = (targetVelocity+aVelocity).mag()/velMag;
00124      random = G4UniformRand();
00125   }
00126   return result;
00127 }
00128 
00129 G4ReactionProduct
00130 G4Nucleus::GetThermalNucleus(G4double targetMass, G4double temp) const
00131 {
00132     G4double currentTemp = temp;
00133     if(currentTemp < 0) currentTemp = theTemp;
00134     G4ReactionProduct theTarget;    
00135     theTarget.SetMass(targetMass*G4Neutron::Neutron()->GetPDGMass());
00136     G4double px, py, pz;
00137     px = GetThermalPz(theTarget.GetMass(), currentTemp);
00138     py = GetThermalPz(theTarget.GetMass(), currentTemp);
00139     pz = GetThermalPz(theTarget.GetMass(), currentTemp);
00140     theTarget.SetMomentum(px, py, pz);
00141     G4double tMom = std::sqrt(px*px+py*py+pz*pz);
00142     G4double tEtot = std::sqrt((tMom+theTarget.GetMass())*
00143                           (tMom+theTarget.GetMass())-
00144                           2.*tMom*theTarget.GetMass());
00145     if(1-tEtot/theTarget.GetMass()>0.001)
00146     {
00147       theTarget.SetTotalEnergy(tEtot);
00148     }
00149     else
00150     {
00151       theTarget.SetKineticEnergy(tMom*tMom/(2.*theTarget.GetMass()));
00152     }    
00153     return theTarget;
00154 }
00155 
00156  
00157 void
00158 G4Nucleus::ChooseParameters(const G4Material* aMaterial)
00159 {
00160   G4double random = G4UniformRand();
00161   G4double sum = aMaterial->GetTotNbOfAtomsPerVolume();
00162   const G4ElementVector* theElementVector = aMaterial->GetElementVector();
00163   G4double running(0);
00164   //  G4Element* element(0);
00165   G4Element* element = (*theElementVector)[aMaterial->GetNumberOfElements()-1];
00166 
00167   for (unsigned int i = 0; i < aMaterial->GetNumberOfElements(); ++i) {
00168     running += aMaterial->GetVecNbOfAtomsPerVolume()[i];
00169     if (running > random*sum) {
00170       element = (*theElementVector)[i];
00171       break;
00172     }
00173   }
00174 
00175   if (element->GetNumberOfIsotopes() > 0) {
00176     G4double randomAbundance = G4UniformRand();
00177     G4double sumAbundance = element->GetRelativeAbundanceVector()[0];
00178     unsigned int iso=0;
00179     while (iso < element->GetNumberOfIsotopes() &&
00180            sumAbundance < randomAbundance) {
00181       ++iso;
00182       sumAbundance += element->GetRelativeAbundanceVector()[iso];
00183     }
00184     theA=element->GetIsotope(iso)->GetN();
00185     theZ=element->GetIsotope(iso)->GetZ();
00186     aEff=theA;
00187     zEff=theZ;
00188   } else {   
00189     aEff = element->GetN();
00190     zEff = element->GetZ();
00191     theZ = G4int(zEff + 0.5);
00192     theA = G4int(aEff + 0.5);   
00193   }
00194 }
00195 
00196 
00197 void
00198 G4Nucleus::SetParameters(G4double A, G4double Z)
00199 {
00200   theZ = G4lrint(Z);
00201   theA = G4lrint(A);   
00202   if (theA<1 || theZ<0 || theZ>theA) {
00203     throw G4HadronicException(__FILE__, __LINE__,
00204             "G4Nucleus::SetParameters called with non-physical parameters");
00205   }
00206   aEff = A;  // atomic weight
00207   zEff = Z;  // atomic number
00208   fIsotope = 0;
00209 }
00210 
00211 void
00212 G4Nucleus::SetParameters(G4int A, const G4int Z )
00213 {
00214   theZ = Z;
00215   theA = A;   
00216   if( theA<1 || theZ<0 || theZ>theA )
00217     {
00218       throw G4HadronicException(__FILE__, __LINE__,
00219                                 "G4Nucleus::SetParameters called with non-physical parameters");
00220     }
00221   aEff = A;  // atomic weight
00222   zEff = Z;  // atomic number
00223   fIsotope = 0;
00224 }
00225 
00226  G4DynamicParticle *
00227   G4Nucleus::ReturnTargetParticle() const
00228   {
00229     // choose a proton or a neutron as the target particle
00230     
00231     G4DynamicParticle *targetParticle = new G4DynamicParticle;
00232     if( G4UniformRand() < zEff/aEff )
00233       targetParticle->SetDefinition( G4Proton::Proton() );
00234     else
00235       targetParticle->SetDefinition( G4Neutron::Neutron() );
00236     return targetParticle;
00237   }
00238  
00239  G4double
00240   G4Nucleus::AtomicMass( const G4double A, const G4double Z ) const
00241   {
00242     // Now returns (atomic mass - electron masses) 
00243     return G4NucleiProperties::GetNuclearMass(A, Z);
00244   }
00245  
00246  G4double
00247   G4Nucleus::AtomicMass( const G4int A, const G4int Z ) const
00248   {
00249     // Now returns (atomic mass - electron masses) 
00250     return G4NucleiProperties::GetNuclearMass(A, Z);
00251   }
00252  
00253  G4double
00254   G4Nucleus::GetThermalPz( const G4double mass, const G4double temp ) const
00255   {
00256     G4double result = G4RandGauss::shoot();
00257     result *= std::sqrt(k_Boltzmann*temp*mass); // Das ist impuls (Pz),
00258                                            // nichtrelativistische rechnung
00259                                            // Maxwell verteilung angenommen
00260     return result;
00261   }
00262  
00263  G4double 
00264   G4Nucleus::EvaporationEffects( G4double kineticEnergy )
00265   {
00266     // derived from original FORTRAN code EXNU by H. Fesefeldt (10-Dec-1986)
00267     //
00268     // Nuclear evaporation as function of atomic number
00269     // and kinetic energy (MeV) of primary particle
00270     //
00271     // returns kinetic energy (MeV)
00272     //
00273     if( aEff < 1.5 )
00274     {
00275       pnBlackTrackEnergy = dtaBlackTrackEnergy = 0.0;
00276       return 0.0;
00277     }
00278     G4double ek = kineticEnergy/GeV;
00279     G4float ekin = std::min( 4.0, std::max( 0.1, ek ) );
00280     const G4float atno = std::min( 120., aEff ); 
00281     const G4float gfa = 2.0*((aEff-1.0)/70.)*std::exp(-(aEff-1.0)/70.);
00282     //
00283     // 0.35 value at 1 GeV
00284     // 0.05 value at 0.1 GeV
00285     //
00286     G4float cfa = std::max( 0.15, 0.35 + ((0.35-0.05)/2.3)*std::log(ekin) );
00287     G4float exnu = 7.716 * cfa * std::exp(-cfa)
00288       * ((atno-1.0)/120.)*std::exp(-(atno-1.0)/120.);
00289     G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin*ekin );
00290     //
00291     // pnBlackTrackEnergy  is the kinetic energy (in GeV) available for
00292     //                     proton/neutron black track particles
00293     // dtaBlackTrackEnergy is the kinetic energy (in GeV) available for
00294     //                     deuteron/triton/alpha black track particles
00295     //
00296     pnBlackTrackEnergy = exnu*fpdiv;
00297     dtaBlackTrackEnergy = exnu*(1.0-fpdiv);
00298     
00299     if( G4int(zEff+0.1) != 82 )
00300     { 
00301       G4double ran1 = -6.0;
00302       G4double ran2 = -6.0;
00303       for( G4int i=0; i<12; ++i )
00304       {
00305         ran1 += G4UniformRand();
00306         ran2 += G4UniformRand();
00307       }
00308       pnBlackTrackEnergy *= 1.0 + ran1*gfa;
00309       dtaBlackTrackEnergy *= 1.0 + ran2*gfa;
00310     }
00311     pnBlackTrackEnergy = std::max( 0.0, pnBlackTrackEnergy );
00312     dtaBlackTrackEnergy = std::max( 0.0, dtaBlackTrackEnergy );
00313     while( pnBlackTrackEnergy+dtaBlackTrackEnergy >= ek )
00314     {
00315       pnBlackTrackEnergy *= 1.0 - 0.5*G4UniformRand();
00316       dtaBlackTrackEnergy *= 1.0 - 0.5*G4UniformRand();
00317     }
00318 //    G4cout << "EvaporationEffects "<<kineticEnergy<<" "
00319 //           <<pnBlackTrackEnergy+dtaBlackTrackEnergy<<endl;
00320     return (pnBlackTrackEnergy+dtaBlackTrackEnergy)*GeV;
00321   }
00322  
00323  G4double G4Nucleus::AnnihilationEvaporationEffects(G4double kineticEnergy, G4double ekOrg)
00324   {
00325     // Nuclear evaporation as a function of atomic number and kinetic 
00326     // energy (MeV) of primary particle.  Modified for annihilation effects. 
00327     //
00328     if( aEff < 1.5 || ekOrg < 0.)
00329     {
00330       pnBlackTrackEnergyfromAnnihilation = 0.0;
00331       dtaBlackTrackEnergyfromAnnihilation = 0.0;
00332       return 0.0;
00333     }
00334     G4double ek = kineticEnergy/GeV;
00335     G4float ekin = std::min( 4.0, std::max( 0.1, ek ) );
00336     const G4float atno = std::min( 120., aEff ); 
00337     const G4float gfa = 2.0*((aEff-1.0)/70.)*std::exp(-(aEff-1.0)/70.);
00338 
00339     G4float cfa = std::max( 0.15, 0.35 + ((0.35-0.05)/2.3)*std::log(ekin) );
00340     G4float exnu = 7.716 * cfa * std::exp(-cfa)
00341       * ((atno-1.0)/120.)*std::exp(-(atno-1.0)/120.);
00342     G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin*ekin );
00343 
00344     pnBlackTrackEnergyfromAnnihilation = exnu*fpdiv;
00345     dtaBlackTrackEnergyfromAnnihilation = exnu*(1.0-fpdiv);
00346     
00347     G4double ran1 = -6.0;
00348     G4double ran2 = -6.0;
00349     for( G4int i=0; i<12; ++i ) {
00350       ran1 += G4UniformRand();
00351       ran2 += G4UniformRand();
00352     }
00353     pnBlackTrackEnergyfromAnnihilation *= 1.0 + ran1*gfa;
00354     dtaBlackTrackEnergyfromAnnihilation *= 1.0 + ran2*gfa;
00355 
00356     pnBlackTrackEnergyfromAnnihilation = std::max( 0.0, pnBlackTrackEnergyfromAnnihilation);
00357     dtaBlackTrackEnergyfromAnnihilation = std::max( 0.0, dtaBlackTrackEnergyfromAnnihilation);
00358     G4double blackSum = pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation;
00359     if (blackSum >= ekOrg/GeV) {
00360       pnBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum;
00361       dtaBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum;
00362     }
00363 
00364     return (pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation)*GeV;
00365   }
00366  
00367  G4double 
00368   G4Nucleus::Cinema( G4double kineticEnergy )
00369   {
00370     // derived from original FORTRAN code CINEMA by H. Fesefeldt (14-Oct-1987)
00371     //
00372     // input: kineticEnergy (MeV)
00373     // returns modified kinetic energy (MeV)
00374     //
00375     static const G4double expxu =  82.;           // upper bound for arg. of exp
00376     static const G4double expxl = -expxu;         // lower bound for arg. of exp
00377     
00378     G4double ek = kineticEnergy/GeV;
00379     G4double ekLog = std::log( ek );
00380     G4double aLog = std::log( aEff );
00381     G4double em = std::min( 1.0, 0.2390 + 0.0408*aLog*aLog );
00382     G4double temp1 = -ek * std::min( 0.15, 0.0019*aLog*aLog*aLog );
00383     G4double temp2 = std::exp( std::max( expxl, std::min( expxu, -(ekLog-em)*(ekLog-em)*2.0 ) ) );
00384     G4double result = 0.0;
00385     if( std::abs( temp1 ) < 1.0 )
00386     {
00387       if( temp2 > 1.0e-10 )result = temp1*temp2;
00388     }
00389     else result = temp1*temp2;
00390     if( result < -ek )result = -ek;
00391     return result*GeV;
00392   }
00393 
00394  //
00395  // methods for class G4Nucleus  ... by Christian Volcker
00396  //
00397 
00398  G4ThreeVector G4Nucleus::GetFermiMomentum()
00399   {
00400     // chv: .. we assume zero temperature!
00401     
00402     // momentum is equally distributed in each phasespace volume dpx, dpy, dpz.
00403     G4double ranflat1=
00404       CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);   
00405     G4double ranflat2=
00406       CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);   
00407     G4double ranflat3=
00408       CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);   
00409     G4double ranmax = (ranflat1>ranflat2? ranflat1: ranflat2);
00410     ranmax = (ranmax>ranflat3? ranmax : ranflat3);
00411     
00412     // Isotropic momentum distribution
00413     G4double costheta = 2.*G4UniformRand() - 1.0;
00414     G4double sintheta = std::sqrt(1.0 - costheta*costheta);
00415     G4double phi = 2.0*pi*G4UniformRand();
00416     
00417     G4double pz=costheta*ranmax;
00418     G4double px=sintheta*std::cos(phi)*ranmax;
00419     G4double py=sintheta*std::sin(phi)*ranmax;
00420     G4ThreeVector p(px,py,pz);
00421     return p;
00422   }
00423  
00424  G4ReactionProductVector* G4Nucleus::Fragmentate()
00425   {
00426     // needs implementation!
00427     return NULL;
00428   }
00429  
00430  void G4Nucleus::AddMomentum(const G4ThreeVector aMomentum)
00431   {
00432     momentum+=(aMomentum);
00433   }
00434  
00435  void G4Nucleus::AddExcitationEnergy( G4double anEnergy )
00436   {
00437     excitationEnergy+=anEnergy;
00438   }
00439 
00440  /* end of file */
00441 

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