TestEm7/src/G4LindhardPartition.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file G4LindhardPartition.cc
0027 /// \brief Implementation of the G4LindhardPartition class
0028
0029 /*
0030  *  \file electromagnetic/TestEm7/src/G4LindhardPartition.cc
0031  *  \brief Implementation of the G4LindhardPartition class
0032  *
0033  *  Created by Marcus Mendenhall on 1/14/08.
0034  *  2008 Vanderbilt University, Nashville, TN, USA.
0035  *
0036  */
0037
0038 //
0039
0040 #include "G4LindhardPartition.hh"
0041
0042 #include "G4Element.hh"
0043 #include "G4Material.hh"
0044 #include "G4PhysicalConstants.hh"
0045 #include "G4SystemOfUnits.hh"
0046
0047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0048 /*
0049 for a first cut, we will compute NIEL from a Lindhard-Robinson partition
0050 based on the most abundant element in the material.
0051
0052 this is from IEEE Trans. Nucl Science Vol. 48 No.1 February 2001 page 162++
0053 Insoo Jun, "Effects of Secondary Particles on the Total Dose..."
0054 and, by reference,
0055 Lindhard, Nielsen, Scharff & Thompson,
0056 "Integral Equations Governing Radiation Efects...",
0057 Mat. Fys. Medd. Dan. Vid. Selsk. vol 33 #10, pp1-42, 1963
0058 and
0059 Robinson, "The dependence of radiation effects on primary recoil energy",
0060 in Proc. Int. Conf. Radiation-Induced Voids in Metal,
0061 Albany, NY 1972 pp. 397-439
0062 def lindhard_robinson(z1, a1, z2, a2, ke):
0063 el=30.724*z1*z2*math.sqrt(z1**0.6667+z2**0.6667)*(a1+a2)/a2
0064 fl=0.0793*z1**0.6667*math.sqrt(z2)*(a1+a2)**1.5/
0065 ((z1**0.6667+z2**0.6667)**0.75*a1**1.5*math.sqrt(a2))
0066 eps=ke*(1.0/el)
0067 return 1.0/(1+fl*(3.4008*eps**0.16667+0.40244*eps**0.75+eps))
0068 */
0069
0070 G4LindhardRobinsonPartition::G4LindhardRobinsonPartition()
0071 {
0072   max_z = 120;
0073   for (size_t i = 1; i < max_z; i++) {
0074     z23[i] = std::pow((G4double)i, 2. / 3.);
0075   }
0076 }
0077
0078 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0079
0080 G4double G4LindhardRobinsonPartition::PartitionNIEL(G4int z1, G4double a1,
0081                                                     const G4Material* material,
0082                                                     G4double energy) const
0083 {
0084   size_t nMatElements = material->GetNumberOfElements();
0085
0086   const G4double* atomDensities = material->GetVecNbOfAtomsPerVolume();
0087   G4double maxdens = 0.0;
0088   size_t maxindex = 0;
0089   for (size_t k = 0; k < nMatElements; k++) {
0090     if (atomDensities[k] > maxdens) {
0091       maxdens = atomDensities[k];
0092       maxindex = k;
0093     }
0094   }
0095   const G4Element* element = material->GetElement(maxindex);
0096
0097   G4int z2 = G4int(element->GetZ());
0098
0099   G4double a2 = element->GetA() / (Avogadro * amu);
0100
0101   G4double zpow = z23[z1] + z23[z2];
0102   G4double asum = a1 + a2;
0103
0104   G4double el = 30.724 * z1 * z2 * std::sqrt(zpow) * asum / a2;
0105   G4double fl = 0.0793 * z23[z1] * std::sqrt(z2 * asum * asum * asum / (a1 * a1 * a1 * a2))
0106                 / std::pow(zpow, 0.75);
0107   G4double eps = (energy / eV) * (1.0 / el);
0108
0109   return 1.0 / (1 + fl * (3.4008 * std::pow(eps, 0.16667) + 0.40244 * std::pow(eps, 0.75) + eps));
0110 }