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 "G4NeutronHPFissionFS.hh"
00034 #include "G4PhysicalConstants.hh"
00035 #include "G4Nucleus.hh"
00036 #include "G4DynamicParticleVector.hh"
00037 #include "G4NeutronHPFissionERelease.hh"
00038 #include "G4ParticleTable.hh"
00039
00040 void G4NeutronHPFissionFS::Init (G4double A, G4double Z, G4int M, G4String & dirName, G4String & aFSType)
00041 {
00042
00043 theFS.Init(A, Z, M, dirName, aFSType);
00044 theFC.Init(A, Z, M, dirName, aFSType);
00045 theSC.Init(A, Z, M, dirName, aFSType);
00046 theTC.Init(A, Z, M, dirName, aFSType);
00047 theLC.Init(A, Z, M, dirName, aFSType);
00048
00049 theFF.Init(A, Z, M, dirName, aFSType);
00050 if ( getenv("G4NEUTRONHP_PRODUCE_FISSION_FRAGMENTS") && theFF.HasFSData() )
00051 {
00052 G4cout << "Activate Fission Fragments Prodcution for the target isotope of "
00053 << "Z = " << (G4int)Z
00054 << ", A = " << (G4int)A
00055
00056 << G4endl;
00057 G4cout << "As the result, delayed neutrons are omitted and they should be taken care by RadioaActiveDecay."
00058 << G4endl;
00059 produceFissionFragments = true;
00060 }
00061 }
00062 G4HadFinalState * G4NeutronHPFissionFS::ApplyYourself(const G4HadProjectile & theTrack)
00063 {
00064
00065
00066 theResult.Clear();
00067 G4double eKinetic = theTrack.GetKineticEnergy();
00068 const G4HadProjectile *incidentParticle = &theTrack;
00069 G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) );
00070 theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
00071 theNeutron.SetKineticEnergy( eKinetic );
00072
00073
00074 G4Nucleus aNucleus;
00075 G4ReactionProduct theTarget;
00076 G4double targetMass = theFS.GetMass();
00077 G4ThreeVector neuVelo = (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum();
00078 theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
00079
00080
00081 theFS.SetNeutron(theNeutron);
00082 theFS.SetTarget(theTarget);
00083 theFC.SetNeutron(theNeutron);
00084 theFC.SetTarget(theTarget);
00085 theSC.SetNeutron(theNeutron);
00086 theSC.SetTarget(theTarget);
00087 theTC.SetNeutron(theNeutron);
00088 theTC.SetTarget(theTarget);
00089 theLC.SetNeutron(theNeutron);
00090 theLC.SetTarget(theTarget);
00091
00092
00093 theFF.SetNeutron(theNeutron);
00094 theFF.SetTarget(theTarget);
00095
00096
00097
00098
00099
00101
00102
00103 theNeutron.Lorentz(theNeutron, -1*theTarget);
00104
00105
00106
00107 G4DynamicParticleVector * thePhotons;
00108 thePhotons = theFS.GetPhotons();
00109
00110
00111
00112 eKinetic = theNeutron.GetKineticEnergy();
00113 G4double xSec[4];
00114 xSec[0] = theFC.GetXsec(eKinetic);
00115 xSec[1] = xSec[0]+theSC.GetXsec(eKinetic);
00116 xSec[2] = xSec[1]+theTC.GetXsec(eKinetic);
00117 xSec[3] = xSec[2]+theLC.GetXsec(eKinetic);
00118 G4int it;
00119 unsigned int i=0;
00120 G4double random = G4UniformRand();
00121 if(xSec[3]==0)
00122 {
00123 it=-1;
00124 }
00125 else
00126 {
00127 for(i=0; i<4; i++)
00128 {
00129 it =i;
00130 if(random<xSec[i]/xSec[3]) break;
00131 }
00132 }
00133
00134
00135
00136
00137
00138 G4int Prompt=0, delayed=0, all=0;
00139 G4DynamicParticleVector * theNeutrons = 0;
00140 switch(it)
00141 {
00142 case 0:
00143 theFS.SampleNeutronMult(all, Prompt, delayed, eKinetic, 0);
00144 if(Prompt==0&&delayed==0) Prompt=all;
00145 theNeutrons = theFC.ApplyYourself(Prompt);
00146
00147 break;
00148 case 1:
00149 theFS.SampleNeutronMult(all, Prompt, delayed, eKinetic, 1);
00150 if(Prompt==0&&delayed==0) Prompt=all;
00151 theNeutrons = theSC.ApplyYourself(Prompt);
00152 break;
00153 case 2:
00154 theFS.SampleNeutronMult(all, Prompt, delayed, eKinetic, 2);
00155 if(Prompt==0&&delayed==0) Prompt=all;
00156 theNeutrons = theTC.ApplyYourself(Prompt);
00157 break;
00158 case 3:
00159 theFS.SampleNeutronMult(all, Prompt, delayed, eKinetic, 3);
00160 if(Prompt==0&&delayed==0) Prompt=all;
00161 theNeutrons = theLC.ApplyYourself(Prompt);
00162 break;
00163 default:
00164 break;
00165 }
00166
00167
00168
00169
00170
00171 if ( produceFissionFragments ) delayed=0;
00172
00173 G4double * theDecayConstants;
00174
00175 if( theNeutrons != 0)
00176 {
00177 theDecayConstants = new G4double[delayed];
00178
00179
00180
00181
00182 for(i=0; i<theNeutrons->size(); i++)
00183 {
00184 theResult.AddSecondary(theNeutrons->operator[](i));
00185 }
00186 delete theNeutrons;
00187
00188 G4DynamicParticleVector * theDelayed = 0;
00189
00190 theDelayed = theFS.ApplyYourself(0, delayed, theDecayConstants);
00191 for(i=0; i<theDelayed->size(); i++)
00192 {
00193 G4double time = -std::log(G4UniformRand())/theDecayConstants[i];
00194 time += theTrack.GetGlobalTime();
00195 theResult.AddSecondary(theDelayed->operator[](i));
00196 theResult.GetSecondary(theResult.GetNumberOfSecondaries()-1)->SetTime(time);
00197 }
00198 delete theDelayed;
00199 }
00200 else
00201 {
00202
00203 theFS.SampleNeutronMult(all, Prompt, delayed, eKinetic, 0);
00204 theDecayConstants = new G4double[delayed];
00205 if(Prompt==0&&delayed==0) Prompt=all;
00206 theNeutrons = theFS.ApplyYourself(Prompt, delayed, theDecayConstants);
00207
00208
00209
00210 G4int i0;
00211 for(i0=0; i0<Prompt; i0++)
00212 {
00213 theResult.AddSecondary(theNeutrons->operator[](i0));
00214 }
00215
00216
00217 for(i0=Prompt; i0<Prompt+delayed; i0++)
00218 {
00219 G4double time = -std::log(G4UniformRand())/theDecayConstants[i0-Prompt];
00220 time += theTrack.GetGlobalTime();
00221 theResult.AddSecondary(theNeutrons->operator[](i0));
00222 theResult.GetSecondary(theResult.GetNumberOfSecondaries()-1)->SetTime(time);
00223 }
00224 delete theNeutrons;
00225 }
00226 delete [] theDecayConstants;
00227
00228 unsigned int nPhotons = 0;
00229 if(thePhotons!=0)
00230 {
00231 nPhotons = thePhotons->size();
00232 for(i=0; i<nPhotons; i++)
00233 {
00234 theResult.AddSecondary(thePhotons->operator[](i));
00235 }
00236 delete thePhotons;
00237 }
00238
00239
00240
00241
00242
00243
00244
00245 G4NeutronHPFissionERelease * theERelease = theFS.GetEnergyRelease();
00246 G4double eDepByFragments = theERelease->GetFragmentKinetic();
00247
00248 if ( !produceFissionFragments ) theResult.SetLocalEnergyDeposit(eDepByFragments);
00249
00250
00251 theResult.SetStatusChange(stopAndKill);
00252
00253
00254
00255
00256 if ( produceFissionFragments )
00257 {
00258 G4int fragA_Z=0;
00259 G4int fragA_A=0;
00260 G4int fragA_M=0;
00261
00262 theFF.GetAFissionFragment(eKinetic,fragA_Z,fragA_A,fragA_M);
00263 G4int fragB_Z=(G4int)theBaseZ-fragA_Z;
00264 G4int fragB_A=(G4int)theBaseA-fragA_A-Prompt;
00265
00266
00267
00268
00269
00270 G4ParticleTable* pt = G4ParticleTable::GetParticleTable();
00271
00272 G4ParticleDefinition* pdA = pt->GetIon( fragA_Z , fragA_A , 0.0 );
00273 G4ParticleDefinition* pdB = pt->GetIon( fragB_Z , fragB_A , 0.0 );
00274
00275
00276 G4double phi = twopi*G4UniformRand();
00277 G4double theta = pi*G4UniformRand();
00278 G4double sinth = std::sin(theta);
00279 G4ThreeVector direction (sinth*std::cos(phi) , sinth*std::sin(phi), std::cos(theta) );
00280
00281
00282 G4double ER = eDepByFragments;
00283 G4double ma = pdA->GetPDGMass();
00284 G4double mb = pdB->GetPDGMass();
00285 G4double EA = ER / ( 1 + ma/mb);
00286 G4double EB = ER - EA;
00287 G4DynamicParticle* dpA = new G4DynamicParticle( pdA , direction , EA);
00288 G4DynamicParticle* dpB = new G4DynamicParticle( pdB , -direction , EB);
00289 theResult.AddSecondary(dpA);
00290 theResult.AddSecondary(dpB);
00291 }
00292
00293
00294 return &theResult;
00295 }