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 #include "G4NeutronIsotopeProduction.hh" 00027 #include "G4SystemOfUnits.hh" 00028 00029 G4NeutronIsotopeProduction::G4NeutronIsotopeProduction() 00030 { 00031 numberOfElements = G4Element::GetNumberOfElements(); 00032 theData = new G4ElementIsoCrossSections<G4NeutronIsoIsoCrossSections>* [numberOfElements]; 00033 for (G4int i = 0; i < numberOfElements; i++) { 00034 theData[i] = new G4ElementIsoCrossSections<G4NeutronIsoIsoCrossSections>; 00035 if((*(G4Element::GetElementTable()))[i]->GetZ() > 9 && 00036 (*(G4Element::GetElementTable()))[i]->GetZ() < 84) { 00037 // Workaround to be fixed in G4NeutronHPNames. 00038 theData[i]->Init((*(G4Element::GetElementTable()))[i]); 00039 } 00040 } 00041 G4cout << "WARNING: G4NeutronIsotopeProduction is deprecated and will be removed with Geant4 version 10" 00042 << G4endl; 00043 } 00044 00045 00046 G4NeutronIsotopeProduction::~G4NeutronIsotopeProduction() 00047 { 00048 for (G4int i = 0; i < numberOfElements; i++) delete theData[i]; 00049 if (theData) delete [] theData; 00050 } 00051 00052 00053 G4IsoResult* 00054 G4NeutronIsotopeProduction::GetIsotope(const G4HadProjectile* aTrack, 00055 const G4Nucleus&) 00056 { 00057 G4IsoResult* result = 0; 00058 // is applicable? 00059 if (aTrack->GetDefinition() != G4Neutron::Neutron()) return result; 00060 G4double incidentKE = aTrack->GetKineticEnergy(); 00061 if (incidentKE > 100*MeV) return result; 00062 00063 // get the isotope 00064 const G4Material* theMaterial = aTrack->GetMaterial(); 00065 G4int nEleInMat = theMaterial->GetNumberOfElements(); 00066 for (G4int check = 0; check < nEleInMat; check++) { 00067 // Workaround to be fixed in G4NeutronHPNames 00068 if (theMaterial->GetElement(check)->GetZ() < 10) return result; 00069 // Workaround to be fixed in G4NeutronHPNames. 00070 if (theMaterial->GetElement(check)->GetZ() > 83) return result; 00071 } 00072 G4int index(0); 00073 G4double* xSec = new G4double[nEleInMat]; 00074 G4double sum = 0; 00075 00076 for (G4int i = 0; i < nEleInMat; i++) { 00077 index = theMaterial->GetElement(i)->GetIndex(); 00078 xSec[i] = theData[index]->GetCrossSection(incidentKE); 00079 sum += xSec[i]; 00080 } 00081 00082 G4double random = G4UniformRand(); 00083 G4double running = 0; 00084 00085 for (G4int i = 0; i < nEleInMat; i++) { 00086 running += xSec[i]; 00087 index = theMaterial->GetElement(i)->GetIndex(); 00088 if (random <= running/sum) break; 00089 } 00090 00091 delete [] xSec; 00092 result = theData[index]->GetProductIsotope(incidentKE); 00093 return result; 00094 } 00095