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
00033 #include "G4NeutronHPJENDLHEData.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "G4LPhysicsFreeVector.hh"
00036 #include "G4ElementTable.hh"
00037 #include "G4NeutronHPData.hh"
00038
00039 G4bool G4NeutronHPJENDLHEData::IsApplicable(const G4DynamicParticle*aP, const G4Element* anE)
00040 {
00041
00042 G4bool result = true;
00043 G4double eKin = aP->GetKineticEnergy();
00044
00045 if ( eKin < 20*MeV || 3*GeV < eKin || aP->GetDefinition()!=G4Neutron::Neutron() )
00046 {
00047 result = false;
00048 }
00049
00050 else if ( !(vElement[ anE->GetIndex() ]) ) result = false;
00051
00052 return result;
00053
00054 }
00055
00056
00057
00058 G4NeutronHPJENDLHEData::G4NeutronHPJENDLHEData()
00059 {
00060 ;
00061 }
00062
00063
00064
00065 G4NeutronHPJENDLHEData::G4NeutronHPJENDLHEData( G4String reaction , G4ParticleDefinition* pd )
00066 :G4VCrossSectionDataSet( "JENDLHE"+reaction+"CrossSection" )
00067 {
00068 reactionName = reaction;
00069 BuildPhysicsTable( *pd );
00070 }
00071
00072
00073
00074 G4NeutronHPJENDLHEData::~G4NeutronHPJENDLHEData()
00075 {
00076 ;
00077
00078 }
00079
00080
00081
00082 void G4NeutronHPJENDLHEData::BuildPhysicsTable( const G4ParticleDefinition& aP )
00083 {
00084
00085
00086
00087 particleName = aP.GetParticleName();
00088
00089 G4String baseName = getenv( "G4NEUTRONHPDATA" );
00090 G4String dirName = baseName+"/JENDL_HE/"+particleName+"/"+reactionName ;
00091 G4String aFSType = "/CrossSection/";
00092 G4NeutronHPNames theNames;
00093
00094 G4String filename;
00095
00096
00097
00098
00099 size_t numberOfElements = G4Element::GetNumberOfElements();
00100
00101
00102
00103
00104 static const G4ElementTable *theElementTable = G4Element::GetElementTable();
00105 vElement.clear();
00106 vElement.resize( numberOfElements );
00107 for ( size_t i = 0; i < numberOfElements; ++i )
00108 {
00109
00110 G4Element* theElement = (*theElementTable)[i];
00111 vElement[i] = false;
00112
00113
00114 G4int nIso = (*theElementTable)[i]->GetNumberOfIsotopes();
00115 G4int Z = static_cast<G4int> ((*theElementTable)[i]->GetZ());
00116 if ( nIso!=0 )
00117 {
00118 G4bool found_at_least_one = false;
00119 for ( G4int i1 = 0; i1 < nIso; i1++ )
00120 {
00121 G4int A = theElement->GetIsotope(i1)->GetN();
00122
00123 if ( isThisNewIsotope( Z , A ) )
00124 {
00125
00126 std::stringstream ss;
00127 ss << dirName << aFSType << Z << "_" << A << "_" << theNames.GetName( Z-1 );
00128 filename = ss.str();
00129 std::fstream file;
00130 file.open ( filename , std::fstream::in );
00131 G4int dummy;
00132 file >> dummy;
00133 if ( file.good() )
00134 {
00135
00136
00137 found_at_least_one = true;
00138
00139
00140 G4PhysicsVector* aPhysVec = readAFile ( &file );
00141
00142
00143
00144 registAPhysicsVector( Z , A , aPhysVec );
00145
00146 }
00147 else
00148 {
00149
00150 }
00151
00152 file.close();
00153
00154 }
00155 else
00156 {
00157 found_at_least_one = TRUE;
00158 }
00159 }
00160
00161 if ( found_at_least_one ) vElement[i] = true;
00162
00163 }
00164 else
00165 {
00166 G4StableIsotopes theStableOnes;
00167 G4int first = theStableOnes.GetFirstIsotope( Z );
00168 G4bool found_at_least_one = FALSE;
00169 for ( G4int i1 = 0; i1 < theStableOnes.GetNumberOfIsotopes( static_cast<G4int>(theElement->GetZ() ) ); i1++)
00170 {
00171 G4int A = theStableOnes.GetIsotopeNucleonCount( first+i1 );
00172
00173 if ( isThisNewIsotope( Z , A ) )
00174 {
00175
00176 std::stringstream ss;
00177 ss << dirName << aFSType << Z << "_" << A << "_" << theNames.GetName( Z-1 );
00178 filename = ss.str();
00179
00180 std::fstream file;
00181 file.open ( filename , std::fstream::in );
00182 G4int dummy;
00183 file >> dummy;
00184 if ( file.good() )
00185 {
00186
00187 found_at_least_one = TRUE;
00188
00189
00190 G4PhysicsVector* aPhysVec = readAFile ( &file );
00191
00192
00193 registAPhysicsVector( Z , A , aPhysVec );
00194
00195 }
00196 else
00197 {
00198
00199 }
00200
00201 file.close();
00202 }
00203 else
00204 {
00205 found_at_least_one = TRUE;
00206 }
00207 }
00208
00209 if ( found_at_least_one ) vElement[i] = true;
00210
00211 }
00212
00213 }
00214
00215 }
00216
00217
00218
00219 void G4NeutronHPJENDLHEData::DumpPhysicsTable(const G4ParticleDefinition& aP)
00220 {
00221 if(&aP!=G4Neutron::Neutron())
00222 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
00223
00224 }
00225
00226
00227
00228 G4double G4NeutronHPJENDLHEData::
00229 GetCrossSection(const G4DynamicParticle* aP, const G4Element*anE, G4double )
00230
00231 {
00232
00233
00234
00235
00236
00237
00238
00239 G4double result = 0;
00240
00241 G4double ek = aP->GetKineticEnergy();
00242
00243 G4int nIso = anE->GetNumberOfIsotopes();
00244 G4int Z = static_cast<G4int> ( anE->GetZ() );
00245 if ( nIso!=0 )
00246 {
00247 for ( G4int i1 = 0; i1 < nIso; i1++ )
00248 {
00249
00250 G4int A = anE->GetIsotope(i1)->GetN();
00251 G4double frac = anE->GetRelativeAbundanceVector()[ i1 ];
00252
00253 result += frac * getXSfromThisIsotope( Z , A , ek );
00254
00255
00256 }
00257 }
00258 else
00259 {
00260
00261 G4StableIsotopes theStableOnes;
00262 G4int first = theStableOnes.GetFirstIsotope( Z );
00263 for ( G4int i1 = 0; i1 < theStableOnes.GetNumberOfIsotopes( static_cast<G4int>(anE->GetZ() ) ); i1++)
00264 {
00265
00266 G4int A = theStableOnes.GetIsotopeNucleonCount( first+i1 );
00267 G4double frac = theStableOnes.GetAbundance( first+i1 )*perCent;
00268
00269 result += frac * getXSfromThisIsotope( Z , A , ek );
00270
00271
00272 }
00273 }
00274 return result;
00275
00276 }
00277
00278
00279
00280 G4PhysicsVector* G4NeutronHPJENDLHEData::readAFile ( std::fstream* file )
00281 {
00282
00283 G4int dummy;
00284 G4int len;
00285 *file >> dummy;
00286 *file >> len;
00287
00288 std::vector< G4double > v_e;
00289 std::vector< G4double > v_xs;
00290
00291 for ( G4int i = 0 ; i < len ; i++ )
00292 {
00293 G4double e;
00294 G4double xs;
00295
00296 *file >> e;
00297 *file >> xs;
00298
00299 v_e.push_back( e*eV );
00300 v_xs.push_back( xs*barn );
00301 }
00302
00303 G4LPhysicsFreeVector* aPhysVec = new G4LPhysicsFreeVector( static_cast< size_t >( len ) , v_e.front() , v_e.back() );
00304
00305 for ( G4int i = 0 ; i < len ; i++ )
00306 {
00307 aPhysVec->PutValues( static_cast< size_t >( i ) , v_e[ i ] , v_xs[ i ] );
00308 }
00309
00310 return aPhysVec;
00311 }
00312
00313
00314
00315 G4bool G4NeutronHPJENDLHEData::isThisInMap( G4int z , G4int a )
00316 {
00317 if ( mIsotope.find ( z ) == mIsotope.end() ) return false;
00318 if ( mIsotope.find ( z ) -> second->find ( a ) == mIsotope.find ( z ) -> second->end() ) return false;
00319 return true;
00320 }
00321
00322
00323
00324 void G4NeutronHPJENDLHEData::registAPhysicsVector( G4int Z , G4int A , G4PhysicsVector* aPhysVec )
00325 {
00326
00327 std::pair< G4int , G4PhysicsVector* > aPair = std::pair < G4int , G4PhysicsVector* > ( A , aPhysVec );
00328
00329 std::map < G4int , std::map< G4int , G4PhysicsVector* >* >::iterator itm;
00330 itm = mIsotope.find ( Z );
00331 if ( itm != mIsotope.end() )
00332 {
00333 itm->second->insert ( aPair );
00334 }
00335 else
00336 {
00337 std::map< G4int , G4PhysicsVector* >* aMap = new std::map< G4int , G4PhysicsVector* >;
00338 aMap->insert ( aPair );
00339 mIsotope.insert( std::pair< G4int , std::map< G4int , G4PhysicsVector* >* > ( Z , aMap ) );
00340 }
00341
00342 }
00343
00344
00345
00346 G4double G4NeutronHPJENDLHEData::getXSfromThisIsotope( G4int Z , G4int A , G4double ek )
00347 {
00348
00349 G4double aXSection = 0.0;
00350 G4bool outOfRange;
00351
00352 G4PhysicsVector* aPhysVec;
00353 if ( mIsotope.find ( Z )->second->find ( A ) != mIsotope.find ( Z )->second->end() )
00354 {
00355
00356 aPhysVec = mIsotope.find ( Z )->second->find ( A )->second;
00357 aXSection = aPhysVec->GetValue( ek , outOfRange );
00358
00359 }
00360 else
00361 {
00362
00363
00364 std::map < G4int , G4PhysicsVector* >::iterator it;
00365 G4int delta0 = 99;
00366 for ( it = mIsotope.find ( Z )->second->begin() ; it != mIsotope.find ( Z )->second->end() ; it++ )
00367 {
00368 G4int delta = std::abs( A - it->first );
00369 if ( delta < delta0 ) delta0 = delta;
00370 }
00371
00372
00373 if ( G4UniformRand() < 0.5 ) delta0 *= -1;
00374 G4int A1 = A + delta0;
00375 if ( mIsotope.find ( Z )->second->find ( A1 ) != mIsotope.find ( Z )->second->end() )
00376 {
00377 aPhysVec = mIsotope.find ( Z )->second->find ( A1 )->second;
00378 }
00379 else
00380 {
00381 A1 = A - delta0;
00382 aPhysVec = mIsotope.find ( Z )->second->find ( A1 )->second;
00383 }
00384
00385 aXSection = aPhysVec->GetValue( ek , outOfRange );
00386
00387 aXSection *= std::pow ( 1.0*A/ A1 , 2.0 / 3.0 );
00388
00389 }
00390
00391 return aXSection;
00392 }