G4IonsKoxCrossSection.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 // 18-Sep-2003 First version is written by T. Koi
00027 // 10-Nov-2003 Bug fix at Cal. ke_per_n and D T. Koi
00028 // 12-Nov-2003 Add energy check at lower side T. Koi
00029 // 26-Dec-2006 Add isotope dependence D. Wright
00030 // 14-Mar-2011 Moved constructor, destructor and virtual methods to source by V.Ivanchenko
00031 // 19-Aug-2011 V.Ivanchenko move to new design and make x-section per element
00032 
00033 #include "G4IonsKoxCrossSection.hh"
00034 #include "G4PhysicalConstants.hh"
00035 #include "G4SystemOfUnits.hh"
00036 #include "G4DynamicParticle.hh"
00037 #include "G4NucleiProperties.hh"
00038 #include "G4HadTmpUtil.hh"
00039 #include "G4NistManager.hh"
00040 
00041 G4IonsKoxCrossSection::G4IonsKoxCrossSection()
00042   : G4VCrossSectionDataSet("IonsKox"), lowerLimit ( 10*MeV ), 
00043     r0 ( 1.1*fermi ), rc ( 1.3*fermi )
00044 {}
00045 
00046 G4IonsKoxCrossSection::~G4IonsKoxCrossSection()
00047 {}
00048 
00049 void
00050 G4IonsKoxCrossSection::CrossSectionDescription(std::ostream& outFile) const
00051 {
00052   outFile << "G4IonsKoxCrossSection calculates the total reaction cross\n"
00053           << "section for nucleus-nucleus scattering using the Kox\n"
00054           << "parameterization.  It is valid for projectiles and targets\n"
00055           << "of all Z, at projectile energies up to 10 GeV/n.  If the\n"
00056           << "projectile energy is less than 10 MeV/n, a zero cross section\n"
00057           << "is returned.\n";
00058 }
00059 
00060 G4bool G4IonsKoxCrossSection::IsElementApplicable(const G4DynamicParticle* aDP, 
00061                                                   G4int, const G4Material*)
00062 {
00063   return (1 <= aDP->GetDefinition()->GetBaryonNumber());
00064 }
00065 
00066 G4double 
00067 G4IonsKoxCrossSection::GetElementCrossSection(
00068   const G4DynamicParticle* aParticle, G4int ZZ, const G4Material*)
00069 {
00070    G4double xsection = 0.0;
00071 
00072    G4int Ap = aParticle->GetDefinition()->GetBaryonNumber();
00073    G4int Zp = G4int(aParticle->GetDefinition()->GetPDGCharge() / eplus + 0.5); 
00074    G4double ke_per_N = aParticle->GetKineticEnergy() / Ap;
00075 
00076    // Apply energy check, if less than lower limit then 0 value is returned
00077    // if (  ke_per_N < lowerLimit ) return xsection;
00078 
00079    G4int At = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(ZZ));
00080    G4int Zt = ZZ;  
00081 
00082    G4double one_third = 1.0 / 3.0;
00083 
00084    G4double cubicrAt = std::pow ( G4double(At) , G4double(one_third) );  
00085    G4double cubicrAp = std::pow ( G4double(Ap) , G4double(one_third) );  
00086 
00087    // rc divide fermi
00088    G4double Bc = Zt * Zp / ( (rc/fermi) * (cubicrAp+cubicrAt) );
00089 
00090    G4double targ_mass = G4NucleiProperties::GetNuclearMass(At, Zt);
00091    G4double proj_mass = aParticle->GetMass();
00092    G4double proj_momentum = aParticle->GetMomentum().mag(); 
00093 
00094    G4double Ecm = calEcm ( proj_mass , targ_mass , proj_momentum ); 
00095    if( Ecm <= Bc) return xsection;
00096 
00097    G4double Rvol = r0 * (  cubicrAp + cubicrAt );
00098 
00099 //   G4double ke_per_N = aParticle->GetKineticEnergy() / Ap;
00100    G4double c = calCeValue ( ke_per_N / MeV  );  
00101 
00102    G4double a = 1.85;
00103    G4double Rsurf = r0 * (a*cubicrAp * cubicrAt/(cubicrAp + cubicrAt) - c); 
00104    G4double D = 5.0 * ( At - 2 * Zt ) * Zp / ( Ap * At );
00105    Rsurf = Rsurf + D * fermi;  // multiply D by fermi 
00106 
00107    G4double Rint = Rvol + Rsurf;
00108    xsection = pi * Rint * Rint * ( 1 - Bc / ( Ecm / MeV ) );
00109   
00110    return xsection; 
00111 }
00112 
00113 G4double
00114 G4IonsKoxCrossSection::calEcm(G4double mp, G4double mt, G4double Plab)
00115 {
00116    G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
00117    G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
00118    G4double Pcm = Plab * mt / Ecm;
00119    G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp;
00120    return KEcm;
00121 }
00122 
00123 
00124 G4double G4IonsKoxCrossSection::calCeValue(const G4double ke)
00125 {
00126    // Calculate c value 
00127    // This value is indepenent from projectile and target particle 
00128    // ke is projectile kinetic energy per nucleon in the Lab system with MeV unit 
00129    // fitting function is made by T. Koi 
00130    // There are no data below 30 MeV/n in Kox et al., 
00131 
00132    G4double Ce; 
00133    G4double log10_ke = std::log10 ( ke );   
00134    if (log10_ke > 1.5) 
00135    {
00136       Ce = - 10.0 / std::pow ( G4double(log10_ke) , G4double(5) ) + 2.0;
00137    }
00138    else
00139    {
00140       Ce = (-10.0/std::pow(G4double(1.5), G4double(5) ) + 2.0) /
00141            std::pow(G4double(1.5), G4double(3)) * std::pow(G4double(log10_ke), G4double(3) );
00142 
00143    }
00144    return Ce;
00145 }
00146 

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