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 #include "G4NeutronHPFFFissionFS.hh"
00031 #include "G4SystemOfUnits.hh"
00032
00033 void G4NeutronHPFFFissionFS::Init (G4double A, G4double Z, G4int M, G4String & dirName, G4String & )
00034 {
00035
00036 G4String aString = "FF";
00037
00038 G4String tString = dirName;
00039 G4bool dbool;
00040 G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, tString, aString , dbool);
00041 G4String filename = aFile.GetName();
00042 theBaseA = aFile.GetA();
00043 theBaseZ = aFile.GetZ();
00044
00045
00046 if ( !dbool || ( Z < 2.5 && ( std::abs(theBaseZ-Z)>0.0001 || std::abs(theBaseA-A)>0.0001) ) )
00047 {
00048 hasAnyData = false;
00049 hasFSData = false;
00050 hasXsec = false;
00051 return;
00052 }
00053 std::ifstream theData(filename, std::ios::in);
00054 G4double dummy;
00055 if ( !theData )
00056 {
00057 theData.close();
00058 hasFSData = false;
00059 hasXsec = false;
00060 hasAnyData = false;
00061 return;
00062 }
00063
00064
00065 hasFSData = true;
00066
00067
00068 while ( theData.good() )
00069 {
00070 G4int iMT, iMF;
00071 G4int imax;
00072
00073
00074 theData >> iMT >> iMF >> dummy;
00075
00076 theData >> imax;
00077
00078
00079 std::map< G4double , std::map< G4int , G4double >* >* mEnergyFSPData = new std::map< G4double , std::map< G4int , G4double >* >;
00080
00081 std::map< G4double , G4int >* mInterporation = new std::map< G4double , G4int >;
00082 for ( G4int i = 0 ; i <= imax ; i++ )
00083 {
00084
00085 G4double YY=0.0;
00086 G4double Ei;
00087 G4int jmax;
00088 G4int ip;
00089
00090 theData >> Ei;
00091
00092 theData >> jmax;
00093
00094 theData >> ip;
00095 mInterporation->insert( std::pair<G4double,G4int>(Ei*eV,ip) );
00096
00097 std::map<G4int,G4double>* mFSPYieldData = new std::map<G4int,G4double>;
00098 for ( G4int j = 0 ; j < jmax ; j++ )
00099 {
00100 G4int FSP;
00101 G4int mFSP;
00102 G4double Y;
00103 theData >> FSP >> mFSP >> Y;
00104 G4int k = FSP*100+mFSP;
00105 YY = YY + Y;
00106
00107 mFSPYieldData->insert( std::pair<G4int,G4double>( k , YY ) );
00108 }
00109 mEnergyFSPData->insert( std::pair<G4double,std::map<G4int,G4double>*>(Ei*eV,mFSPYieldData) );
00110 }
00111
00112 FissionProductYieldData.insert( std::pair< G4int , std::map< G4double , std::map< G4int , G4double >* >* > (iMT,mEnergyFSPData));
00113 mMTInterpolation.insert( std::pair<G4int,std::map<G4double,G4int>*> (iMT,mInterporation) );
00114 }
00115 theData.close();
00116 }
00117
00118 G4DynamicParticleVector * G4NeutronHPFFFissionFS::ApplyYourself(G4int nNeutrons)
00119 {
00120 G4DynamicParticleVector * aResult;
00121
00122 aResult = G4NeutronHPFissionBaseFS::ApplyYourself(nNeutrons);
00123 return aResult;
00124 }
00125
00126 void G4NeutronHPFFFissionFS::GetAFissionFragment( G4double energy , G4int& fragZ , G4int& fragA , G4int& fragM )
00127 {
00128
00129
00130 G4double rand =G4UniformRand();
00131
00132
00133 std::map< G4double , std::map< G4int , G4double >* >* mEnergyFSPData = FissionProductYieldData.find( 454 )->second;
00134
00135
00136
00137
00138 G4double key_energy = DBL_MAX;
00139 if ( mEnergyFSPData->size() == 1 )
00140 {
00141 key_energy = mEnergyFSPData->begin()->first;
00142 }
00143 else
00144 {
00145
00146 G4double Dmin=DBL_MAX;
00147 G4int i = 0;
00148 for ( std::map< G4double , std::map< G4int , G4double >* >::iterator it = mEnergyFSPData->begin() ;
00149 it != mEnergyFSPData->end() ; it++ )
00150 {
00151 G4double e = (it->first);
00152 G4double d = std::fabs ( energy - e );
00153 if ( d < Dmin )
00154 {
00155 Dmin = d;
00156 key_energy = e;
00157 }
00158 i++;
00159 }
00160 }
00161
00162 std::map<G4int,G4double>* mFSPYieldData = (*mEnergyFSPData)[key_energy];
00163
00164 G4int ifrag=0;
00165 G4double ceilling = mFSPYieldData->rbegin()->second;
00166 for ( std::map<G4int,G4double>::iterator it = mFSPYieldData->begin() ; it != mFSPYieldData->end() ; it++ )
00167 {
00168
00169 if ( rand <= it->second/ceilling )
00170 {
00171
00172 ifrag = it->first;
00173 break;
00174 }
00175 }
00176
00177 fragZ = ifrag/100000;
00178 fragA = (ifrag%100000)/100;
00179 fragM = (ifrag%100);
00180
00181
00182 }