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
00027
00028
00029
00030
00031
00032 #include "G4NeutronHPElementData.hh"
00033 #include "G4SystemOfUnits.hh"
00034
00035 G4NeutronHPElementData::G4NeutronHPElementData()
00036 {
00037 precision = 0.02;
00038 theFissionData = new G4NeutronHPVector;
00039 theCaptureData = new G4NeutronHPVector;
00040 theElasticData = new G4NeutronHPVector;
00041 theInelasticData = new G4NeutronHPVector;
00042 theIsotopeWiseData = 0;
00043 }
00044
00045 G4NeutronHPElementData::~G4NeutronHPElementData()
00046 {
00047 delete theFissionData;
00048 delete theCaptureData;
00049 delete theElasticData;
00050 delete theInelasticData;
00051 delete [] theIsotopeWiseData;
00052 }
00053
00054 void G4NeutronHPElementData::Init(G4Element * theElement)
00055 {
00056 G4int count = theElement->GetNumberOfIsotopes();
00057 if(count == 0) count +=
00058 theStableOnes.GetNumberOfIsotopes(static_cast<G4int>(theElement->GetZ()));
00059 theIsotopeWiseData = new G4NeutronHPIsoData[count];
00060
00061 count = 0;
00062 G4int nIso = theElement->GetNumberOfIsotopes();
00063 G4int Z = static_cast<G4int> (theElement->GetZ());
00064
00065 if(nIso!=0)
00066 {
00067 for (G4int i1=0; i1<nIso; i1++)
00068 {
00069
00070 G4int A = theElement->GetIsotope(i1)->GetN();
00071 G4int M = theElement->GetIsotope(i1)->Getm();
00072 G4double frac = theElement->GetRelativeAbundanceVector()[i1]/perCent;
00073
00074 UpdateData(A, Z, M, count++, frac);
00075 }
00076 }else{
00077
00078 G4int first = theStableOnes.GetFirstIsotope(Z);
00079
00080 for(G4int i1=0;
00081 i1<theStableOnes.GetNumberOfIsotopes(static_cast<G4int>(theElement->GetZ()) );
00082 i1++)
00083 {
00084
00085 G4int A = theStableOnes.GetIsotopeNucleonCount(first+i1);
00086 G4double frac = theStableOnes.GetAbundance(first+i1);
00087
00088 UpdateData(A, Z, count++, frac);
00089 }
00090 }
00091 theElasticData->ThinOut(precision);
00092 theInelasticData->ThinOut(precision);
00093 theCaptureData->ThinOut(precision);
00094 theFissionData->ThinOut(precision);
00095 }
00096
00097
00098 void G4NeutronHPElementData::UpdateData(G4int A, G4int Z, G4int M, G4int index, G4double abundance)
00099 {
00100
00101
00102
00103 theIsotopeWiseData[index].Init(A, Z, M, abundance);
00104
00105
00106 theBuffer = theIsotopeWiseData[index].MakeElasticData();
00107
00108
00109 Harmonise(theElasticData, theBuffer);
00110
00111
00112 delete theBuffer;
00113
00114 theBuffer = theIsotopeWiseData[index].MakeInelasticData();
00115
00116
00117 Harmonise(theInelasticData, theBuffer);
00118
00119
00120 delete theBuffer;
00121
00122 theBuffer = theIsotopeWiseData[index].MakeCaptureData();
00123
00124
00125 Harmonise(theCaptureData, theBuffer);
00126
00127
00128 delete theBuffer;
00129
00130 theBuffer = theIsotopeWiseData[index].MakeFissionData();
00131
00132
00133 Harmonise(theFissionData, theBuffer);
00134
00135
00136 delete theBuffer;
00137
00138
00139 }
00140
00141 void G4NeutronHPElementData::Harmonise(G4NeutronHPVector *& theStore, G4NeutronHPVector * theNew)
00142 {
00143 if(theNew == 0) { return; }
00144 G4int s_tmp = 0, n=0, m_tmp=0;
00145 G4NeutronHPVector * theMerge = new G4NeutronHPVector(theStore->GetVectorLength());
00146
00147 while ( theStore->GetEnergy(s_tmp)<theNew->GetEnergy(0)&&s_tmp<theStore->GetVectorLength() )
00148 {
00149 theMerge->SetData(m_tmp++, theStore->GetEnergy(s_tmp), theStore->GetXsec(s_tmp));
00150 s_tmp++;
00151 }
00152 G4NeutronHPVector *active = theStore;
00153 G4NeutronHPVector * passive = theNew;
00154 G4NeutronHPVector * tmp;
00155 G4int a = s_tmp, p = n, t;
00156
00157 while (a<active->GetVectorLength()&&p<passive->GetVectorLength())
00158 {
00159 if(active->GetEnergy(a) <= passive->GetEnergy(p))
00160 {
00161 theMerge->SetData(m_tmp, active->GetEnergy(a), active->GetXsec(a));
00162 G4double x = theMerge->GetEnergy(m_tmp);
00163 G4double y = std::max(0., passive->GetXsec(x));
00164 theMerge->SetData(m_tmp, x, theMerge->GetXsec(m_tmp)+y);
00165 m_tmp++;
00166 a++;
00167 } else {
00168
00169 tmp = active; t=a;
00170 active = passive; a=p;
00171 passive = tmp; p=t;
00172 }
00173 }
00174
00175 while (a!=active->GetVectorLength())
00176 {
00177 theMerge->SetData(m_tmp++, active->GetEnergy(a), active->GetXsec(a));
00178 a++;
00179 }
00180
00181 while (p!=passive->GetVectorLength())
00182 {
00183
00184
00185 G4double x = passive->GetEnergy(p);
00186 G4double y = std::max(0., active->GetXsec(x));
00187 theMerge->SetData(m_tmp++, x, passive->GetXsec(p)+y);
00188 p++;
00189 }
00190
00191 delete theStore;
00192 theStore = theMerge;
00193
00194 }
00195
00196 G4NeutronHPVector * G4NeutronHPElementData::MakePhysicsVector(G4Element * theElement,
00197 G4ParticleDefinition * theP,
00198 G4NeutronHPFissionData* theSet)
00199 {
00200 if(theP != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
00201 Init ( theElement );
00202 return GetData(theSet);
00203 }
00204 G4NeutronHPVector * G4NeutronHPElementData::MakePhysicsVector(G4Element * theElement,
00205 G4ParticleDefinition * theP,
00206 G4NeutronHPCaptureData * theSet)
00207 {
00208 if(theP != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
00209 Init ( theElement );
00210 return GetData(theSet);
00211 }
00212 G4NeutronHPVector * G4NeutronHPElementData::MakePhysicsVector(G4Element * theElement,
00213 G4ParticleDefinition * theP,
00214 G4NeutronHPElasticData * theSet)
00215 {
00216 if(theP != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
00217 Init ( theElement );
00218 return GetData(theSet);
00219 }
00220 G4NeutronHPVector * G4NeutronHPElementData::MakePhysicsVector(G4Element * theElement,
00221 G4ParticleDefinition * theP,
00222 G4NeutronHPInelasticData * theSet)
00223 {
00224 if(theP != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
00225 Init ( theElement );
00226 return GetData(theSet);
00227 }