00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "G4NeutronIsoIsoCrossSections.hh"
00027 #include "G4SystemOfUnits.hh"
00028 #include "G4NeutronHPDataUsed.hh"
00029 #include "G4NeutronInelasticCrossSection.hh"
00030
00031 G4NeutronIsoIsoCrossSections::
00032 G4NeutronIsoIsoCrossSections()
00033 : theCrossSection(), theNames()
00034 {
00035 theProductionData = NULL;
00036 hasData = false;
00037 theNumberOfProducts = 0;
00038 theZ = 0;
00039 theA = 0;
00040 G4cout << "WARNING: G4NeutronIsoIsoCrossSections is deprecated and will be removed with Geant4 version 10"
00041 << G4endl;
00042 }
00043
00044 G4NeutronIsoIsoCrossSections::
00045 ~G4NeutronIsoIsoCrossSections()
00046 {
00047 for(G4int i=0; i<theNumberOfProducts; i++)
00048 {
00049 delete theProductionData[i];
00050 }
00051 delete [] theProductionData;
00052 }
00053
00054 void G4NeutronIsoIsoCrossSections::
00055 Init(G4int A, G4int Z, G4double frac)
00056 {
00057 frac = frac/100.;
00058
00059
00060
00061 theZ = Z;
00062 theA = A;
00063 theNames.SetMaxOffSet(5);
00064 G4NeutronHPDataUsed dataUsed;
00065 G4String rest = "/CrossSection/";
00066 G4String base = getenv("G4NEUTRONHPDATA");
00067 G4String base1 = base + "/Inelastic/";
00068 G4bool hasInelasticData = false;
00069 dataUsed = theNames.GetName(A, Z, base1, rest, hasInelasticData);
00070 G4String aName = dataUsed.GetName();
00071 G4NeutronHPVector inelasticData;
00072 G4double dummy;
00073 G4int total;
00074 if(hasInelasticData)
00075 {
00076 std::ifstream aDataSet(aName, std::ios::in);
00077 aDataSet >> dummy >> dummy >> total;
00078 inelasticData.Init(aDataSet, total, eV);
00079 }
00080 base1 = base + "/Capture/";
00081 G4bool hasCaptureData = false;
00082 dataUsed = theNames.GetName(A, Z, base1, rest, hasCaptureData);
00083 aName = dataUsed.GetName();
00084 G4NeutronHPVector captureData;
00085 if(hasCaptureData)
00086 {
00087 std::ifstream aDataSet(aName, std::ios::in);
00088 aDataSet >> dummy >> dummy >> total;
00089 captureData.Init(aDataSet, total, eV);
00090 }
00091 base1 = base + "/Elastic/";
00092 G4bool hasElasticData = false;
00093 dataUsed = theNames.GetName(A, Z, base1, rest, hasElasticData);
00094 aName = dataUsed.GetName();
00095 G4NeutronHPVector elasticData;
00096 if(hasElasticData)
00097 {
00098 std::ifstream aDataSet(aName, std::ios::in);
00099 aDataSet >> dummy >> dummy >> total;
00100 elasticData.Init(aDataSet, total, eV);
00101 }
00102 base1 = base + "/Fission/";
00103 G4bool hasFissionData = false;
00104 if(Z>=91)
00105 {
00106 dataUsed = theNames.GetName(A, Z, base1, rest, hasFissionData);
00107 aName = dataUsed.GetName();
00108 }
00109 G4NeutronHPVector fissionData;
00110 if(hasFissionData)
00111 {
00112 std::ifstream aDataSet(aName, std::ios::in);
00113 aDataSet >> dummy >> dummy >> total;
00114 fissionData.Init(aDataSet, total, eV);
00115 }
00116 hasData = hasFissionData||hasInelasticData||hasElasticData||hasCaptureData;
00117 G4NeutronHPVector merged, merged1;
00118 if(hasData)
00119 {
00120 if(hasFissionData&&hasInelasticData)
00121 {
00122 merged = fissionData + inelasticData;
00123 }
00124 else if(hasFissionData)
00125 {
00126 merged = fissionData;
00127 }
00128 else if(hasInelasticData)
00129 {
00130 merged = inelasticData;
00131 }
00132
00133 if(hasElasticData&&hasCaptureData)
00134 {
00135 merged1=elasticData + captureData;
00136 }
00137 else if(hasCaptureData)
00138 {
00139 merged1 = captureData;
00140 }
00141 else if(hasElasticData)
00142 {
00143 merged1 = elasticData;
00144 }
00145
00146 if((hasElasticData||hasCaptureData)&&(hasFissionData||hasInelasticData))
00147 {
00148 theCrossSection = merged + merged1;
00149 }
00150 else if(hasElasticData||hasCaptureData)
00151 {
00152 theCrossSection = merged1;
00153 }
00154 else if(hasFissionData||hasInelasticData)
00155 {
00156 theCrossSection = merged;
00157 }
00158 theCrossSection.Times(frac);
00159 }
00160
00161
00162 theNames.SetMaxOffSet(1);
00163 rest = "/CrossSection/";
00164 base1 = base + "/IsotopeProduction/";
00165 G4bool hasIsotopeProductionData;
00166 dataUsed = theNames.GetName(A, Z, base1, rest, hasIsotopeProductionData);
00167 aName = dataUsed.GetName();
00168 if(hasIsotopeProductionData)
00169 {
00170 std::ifstream aDataSet(aName, std::ios::in);
00171 aDataSet>>theNumberOfProducts;
00172 theProductionData = new G4IsoProdCrossSections * [theNumberOfProducts];
00173 for(G4int i=0; i<theNumberOfProducts; i++)
00174 {
00175 G4String dName;
00176 aDataSet >> dName;
00177 aDataSet >> dummy >> dummy;
00178 theProductionData[i] = new G4IsoProdCrossSections(dName);
00179 theProductionData[i]->Init(aDataSet);
00180 }
00181 }
00182 else
00183 {
00184 hasData = false;
00185 }
00186 G4NeutronInelasticCrossSection theParametrization;
00187 G4double lastEnergy = theCrossSection.GetX(theCrossSection.GetVectorLength()-1);
00188 G4double lastValue = theCrossSection.GetY(theCrossSection.GetVectorLength()-1);
00189 G4double norm = theParametrization.GetCrossSection(lastEnergy, Z, A);
00190 G4double increment = 1*MeV;
00191 while(lastEnergy+increment<101*MeV)
00192 {
00193 G4double currentEnergy = lastEnergy+increment;
00194 G4double value = theParametrization.GetCrossSection(currentEnergy, Z, A)*(lastValue/norm);
00195 G4int position = theCrossSection.GetVectorLength();
00196 theCrossSection.SetData(position, currentEnergy, value);
00197 increment+=1*MeV;
00198 }
00199 }
00200
00201 G4double G4NeutronIsoIsoCrossSections::
00202 GetCrossSection(G4double anEnergy)
00203 {
00204 G4double result;
00205 result = theCrossSection.GetY(anEnergy);
00206 return result;
00207 }
00208
00209 G4String G4NeutronIsoIsoCrossSections::
00210 GetProductIsotope(G4double anEnergy)
00211 {
00212 G4String result = "UNCHANGED";
00213
00214 G4double totalXSec = theCrossSection.GetY(anEnergy);
00215 if(totalXSec==0) return result;
00216
00217
00218 G4double * xSec = new G4double[theNumberOfProducts];
00219 G4double sum = 0;
00220 {
00221 for(G4int i=0; i<theNumberOfProducts; i++)
00222 {
00223 xSec[i] = theProductionData[i]->GetProductionCrossSection(anEnergy);
00224 sum += xSec[i];
00225 }
00226 }
00227 G4double isoChangeXsec = sum;
00228 G4double rand = G4UniformRand();
00229 if(rand > isoChangeXsec/totalXSec)
00230 {
00231 delete [] xSec;
00232 return result;
00233 }
00234 G4double random = G4UniformRand();
00235 G4double running = 0;
00236 G4int index(0);
00237 {
00238 for(G4int i=0; i<theNumberOfProducts; i++)
00239 {
00240 running += xSec[i];
00241 index = i;
00242 if(random<=running/sum) break;
00243 }
00244 }
00245 delete [] xSec;
00246 result = theProductionData[index]->GetProductIsotope();
00247
00248 return result;
00249 }