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 "G4NeutronHPFSFissionFS.hh"
00031 #include "G4ReactionProduct.hh"
00032 #include "G4Nucleus.hh"
00033 #include "G4Proton.hh"
00034 #include "G4Deuteron.hh"
00035 #include "G4Triton.hh"
00036 #include "G4Alpha.hh"
00037 #include "G4ThreeVector.hh"
00038 #include "G4Poisson.hh"
00039 #include "G4LorentzVector.hh"
00040 #include "G4NeutronHPDataUsed.hh"
00041
00042 void G4NeutronHPFSFissionFS::Init (G4double A, G4double Z, G4int M, G4String & dirName, G4String & )
00043 {
00044 G4String tString = "/FS/";
00045 G4bool dbool;
00046 G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, dirName, tString, dbool);
00047 G4String filename = aFile.GetName();
00048 SetAZMs( A, Z, M, aFile );
00049 if(!dbool)
00050 {
00051 hasAnyData = false;
00052 hasFSData = false;
00053 hasXsec = false;
00054 return;
00055 }
00056 std::ifstream theData(filename, std::ios::in);
00057
00058
00059 G4int infoType, dataType;
00060 hasFSData = false;
00061 while (theData >> infoType)
00062 {
00063 hasFSData = true;
00064 theData >> dataType;
00065 switch(infoType)
00066 {
00067 case 1:
00068 if(dataType==4) theNeutronAngularDis.Init(theData);
00069 if(dataType==5) thePromptNeutronEnDis.Init(theData);
00070 if(dataType==12) theFinalStatePhotons.InitMean(theData);
00071 if(dataType==14) theFinalStatePhotons.InitAngular(theData);
00072 if(dataType==15) theFinalStatePhotons.InitEnergies(theData);
00073 break;
00074 case 2:
00075 if(dataType==1) theFinalStateNeutrons.InitMean(theData);
00076 break;
00077 case 3:
00078 if(dataType==1) theFinalStateNeutrons.InitDelayed(theData);
00079 if(dataType==5) theDelayedNeutronEnDis.Init(theData);
00080 break;
00081 case 4:
00082 if(dataType==1) theFinalStateNeutrons.InitPrompt(theData);
00083 break;
00084 case 5:
00085 if(dataType==1) theEnergyRelease.Init(theData);
00086 break;
00087 default:
00088 G4cout << "G4NeutronHPFSFissionFS::Init: unknown data type"<<dataType<<G4endl;
00089 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPFSFissionFS::Init: unknown data type");
00090 break;
00091 }
00092 }
00093 targetMass = theFinalStateNeutrons.GetTargetMass();
00094 theData.close();
00095 }
00096
00097
00098 G4DynamicParticleVector * G4NeutronHPFSFissionFS::ApplyYourself(G4int nPrompt,
00099 G4int nDelayed, G4double * theDecayConst)
00100 {
00101 G4int i;
00102 G4DynamicParticleVector * aResult = new G4DynamicParticleVector;
00103 G4ReactionProduct boosted;
00104 boosted.Lorentz(theNeutron, theTarget);
00105 G4double eKinetic = boosted.GetKineticEnergy();
00106
00107
00108 G4ReactionProduct * theNeutrons = new G4ReactionProduct[nPrompt+nDelayed];
00109 for(i=0; i<nPrompt+nDelayed; i++)
00110 {
00111 theNeutrons[i].SetDefinition(G4Neutron::Neutron());
00112 }
00113
00114
00115 G4int it, dummy;
00116 G4double tempE;
00117 for(i=0; i<nPrompt; i++)
00118 {
00119 tempE = thePromptNeutronEnDis.Sample(eKinetic, dummy);
00120 theNeutrons[i].SetKineticEnergy(tempE);
00121 }
00122 for(i=nPrompt; i<nPrompt+nDelayed; i++)
00123 {
00124 theNeutrons[i].SetKineticEnergy(theDelayedNeutronEnDis.Sample(eKinetic, it));
00125 if(it==0) theNeutrons[i].SetKineticEnergy(thePromptNeutronEnDis.Sample(eKinetic, dummy));
00126 theDecayConst[i-nPrompt] = theFinalStateNeutrons.GetDecayConstant(it);
00127 }
00128
00129
00130 for(i=0; i<nPrompt+nDelayed; i++)
00131 {
00132 theNeutronAngularDis.SampleAndUpdate(theNeutrons[i]);
00133 }
00134
00135
00136 for(i=0; i<nPrompt+nDelayed; i++)
00137 {
00138 G4DynamicParticle * dp = new G4DynamicParticle;
00139 dp->SetDefinition(theNeutrons[i].GetDefinition());
00140 dp->SetMomentum(theNeutrons[i].GetMomentum());
00141 aResult->push_back(dp);
00142 }
00143 delete [] theNeutrons;
00144
00145 return aResult;
00146 }
00147
00148 void G4NeutronHPFSFissionFS::SampleNeutronMult(G4int&all, G4int&Prompt, G4int&delayed, G4double eKinetic, G4int off)
00149 {
00150 G4double promptNeutronMulti = 0;
00151 promptNeutronMulti = theFinalStateNeutrons.GetPrompt(eKinetic);
00152 G4double delayedNeutronMulti = 0;
00153 delayedNeutronMulti = theFinalStateNeutrons.GetDelayed(eKinetic);
00154
00155 if(delayedNeutronMulti==0&&promptNeutronMulti==0)
00156 {
00157 Prompt = 0;
00158 delayed = 0;
00159 G4double totalNeutronMulti = theFinalStateNeutrons.GetMean(eKinetic);
00160 all = G4Poisson(totalNeutronMulti-off);
00161 all += off;
00162 }
00163 else
00164 {
00165 Prompt = G4Poisson(promptNeutronMulti-off);
00166 Prompt += off;
00167 delayed = G4Poisson(delayedNeutronMulti);
00168 all = Prompt+delayed;
00169 }
00170 }
00171
00172 G4DynamicParticleVector * G4NeutronHPFSFissionFS::GetPhotons()
00173 {
00174
00175 G4ReactionProductVector * temp;
00176 G4ReactionProduct boosted;
00177
00178 boosted.Lorentz(theNeutron, theTarget);
00179 G4double anEnergy = boosted.GetKineticEnergy();
00180 temp = theFinalStatePhotons.GetPhotons(anEnergy);
00181 if(temp == 0) { return 0; }
00182
00183
00184 unsigned int i;
00185 G4DynamicParticleVector * result = new G4DynamicParticleVector;
00186 for(i=0; i<temp->size(); i++)
00187 {
00188
00189 temp->operator[](i)->Lorentz(*(temp->operator[](i)), -1.*theTarget);
00190 G4DynamicParticle * theOne = new G4DynamicParticle;
00191 theOne->SetDefinition(temp->operator[](i)->GetDefinition());
00192 theOne->SetMomentum(temp->operator[](i)->GetMomentum());
00193 result->push_back(theOne);
00194 delete temp->operator[](i);
00195 }
00196 delete temp;
00197 return result;
00198 }