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
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063 #include "G4FissionLibrary.hh"
00064 #include "G4SystemOfUnits.hh"
00065
00066 G4FissionLibrary::G4FissionLibrary()
00067 : G4NeutronHPFinalState(), theIsotope(0), targetMass(0.0)
00068 {
00069 hasXsec = false;
00070 }
00071
00072 G4FissionLibrary::~G4FissionLibrary()
00073 {}
00074
00075 G4NeutronHPFinalState * G4FissionLibrary::New()
00076 {
00077 G4FissionLibrary * theNew = new G4FissionLibrary;
00078 return theNew;
00079 }
00080
00081
00082 void G4FissionLibrary::Init (G4double A, G4double Z, G4int M, G4String & dirName, G4String &)
00083 {
00084 G4String tString = "/FS/";
00085 G4bool dbool;
00086 theIsotope = static_cast<G4int>(1000*Z+A);
00087 G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, dirName, tString, dbool);
00088 G4String filename = aFile.GetName();
00089
00090 if(!dbool)
00091 {
00092 hasAnyData = false;
00093 hasFSData = false;
00094 hasXsec = false;
00095 return;
00096 }
00097 std::ifstream theData(filename, std::ios::in);
00098
00099
00100 G4int infoType, dataType;
00101 hasFSData = false;
00102 while (theData >> infoType)
00103 {
00104 hasFSData = true;
00105 theData >> dataType;
00106 switch(infoType)
00107 {
00108 case 1:
00109 if(dataType==4) theNeutronAngularDis.Init(theData);
00110 if(dataType==5) thePromptNeutronEnDis.Init(theData);
00111 if(dataType==12) theFinalStatePhotons.InitMean(theData);
00112 if(dataType==14) theFinalStatePhotons.InitAngular(theData);
00113 if(dataType==15) theFinalStatePhotons.InitEnergies(theData);
00114 break;
00115 case 2:
00116 if(dataType==1) theFinalStateNeutrons.InitMean(theData);
00117 break;
00118 case 3:
00119 if(dataType==1) theFinalStateNeutrons.InitDelayed(theData);
00120 if(dataType==5) theDelayedNeutronEnDis.Init(theData);
00121 break;
00122 case 4:
00123 if(dataType==1) theFinalStateNeutrons.InitPrompt(theData);
00124 break;
00125 case 5:
00126 if(dataType==1) theEnergyRelease.Init(theData);
00127 break;
00128 default:
00129 G4cout << "G4FissionLibrary::Init: unknown data type"<<dataType<<G4endl;
00130 throw G4HadronicException(__FILE__, __LINE__, "G4FissionLibrary::Init: unknown data type");
00131 break;
00132 }
00133 }
00134 targetMass = theFinalStateNeutrons.GetTargetMass();
00135 theData.close();
00136 }
00137
00138 G4HadFinalState * G4FissionLibrary::ApplyYourself(const G4HadProjectile & theTrack)
00139 {
00140 theResult.Clear();
00141
00142
00143 G4double eKinetic = theTrack.GetKineticEnergy();
00144 const G4HadProjectile *incidentParticle = &theTrack;
00145 G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) );
00146 theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
00147 theNeutron.SetKineticEnergy( eKinetic );
00148
00149
00150 G4Nucleus aNucleus;
00151 G4ReactionProduct theTarget;
00152 G4ThreeVector neuVelo = (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum();
00153 theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
00154
00155
00156 theNeutronAngularDis.SetNeutron(theNeutron);
00157 theNeutronAngularDis.SetTarget(theTarget);
00158
00159
00160 theNeutron.Lorentz(theNeutron, -1*theTarget);
00161
00162 eKinetic = theNeutron.GetKineticEnergy();
00163
00164
00165
00166
00167
00168 G4int nPrompt=0, gPrompt=0;
00169 SampleMult(theTrack, &nPrompt, &gPrompt, eKinetic);
00170
00171
00172 G4double momentum;
00173 for(G4int i=0; i<nPrompt; i++)
00174 {
00175 G4DynamicParticle * it = new G4DynamicParticle;
00176 it->SetDefinition(G4Neutron::Neutron());
00177 it->SetKineticEnergy(getneng_(&i)*MeV);
00178 momentum = it->GetTotalMomentum();
00179 G4ThreeVector temp(momentum*getndircosu_(&i),
00180 momentum*getndircosv_(&i),
00181 momentum*getndircosw_(&i));
00182 it->SetMomentum( temp );
00183
00184 theResult.AddSecondary(it);
00185
00186 }
00187
00188
00189 for(G4int i=0; i<gPrompt; i++)
00190 {
00191 G4ReactionProduct * thePhoton = new G4ReactionProduct;
00192 thePhoton->SetDefinition(G4Gamma::Gamma());
00193 thePhoton->SetKineticEnergy(getpeng_(&i)*MeV);
00194 momentum = thePhoton->GetTotalMomentum();
00195 G4ThreeVector temp(momentum*getpdircosu_(&i),
00196 momentum*getpdircosv_(&i),
00197 momentum*getpdircosw_(&i));
00198 thePhoton->SetMomentum( temp );
00199 thePhoton->Lorentz(*thePhoton, -1.*theTarget);
00200
00201 G4DynamicParticle * it = new G4DynamicParticle;
00202 it->SetDefinition(thePhoton->GetDefinition());
00203 it->SetMomentum(thePhoton->GetMomentum());
00204
00205
00206 theResult.AddSecondary(it);
00207 delete thePhoton;
00208 }
00209
00210
00211
00212
00213
00214 G4double eDepByFragments = theEnergyRelease.GetFragmentKinetic();
00215 theResult.SetLocalEnergyDeposit(eDepByFragments);
00216
00217
00218 theResult.SetStatusChange(stopAndKill);
00219 return &theResult;
00220 }
00221
00222 void G4FissionLibrary::SampleMult(const G4HadProjectile & theTrack, G4int* nPrompt,
00223 G4int* gPrompt, G4double eKinetic)
00224 {
00225 G4double promptNeutronMulti = 0;
00226 promptNeutronMulti = theFinalStateNeutrons.GetPrompt(eKinetic);
00227 G4double delayedNeutronMulti = 0;
00228 delayedNeutronMulti = theFinalStateNeutrons.GetDelayed(eKinetic);
00229
00230 G4double time = theTrack.GetGlobalTime()/second;
00231 if(delayedNeutronMulti==0&&promptNeutronMulti==0) {
00232
00233
00234
00235 G4double totalNeutronMulti = theFinalStateNeutrons.GetMean(eKinetic);
00236 genfissevt_(&theIsotope, &time, &totalNeutronMulti, &eKinetic);
00237 } else {
00238
00239 genfissevt_(&theIsotope, &time, &promptNeutronMulti, &eKinetic);
00240 }
00241 *nPrompt = getnnu_();
00242 if (*nPrompt == -1) *nPrompt = 0;
00243 *gPrompt = getpnu_();
00244 if (*gPrompt == -1) *gPrompt = 0;
00245 }
00246