#include <G4NeutronHPCaptureFS.hh>
Inheritance diagram for G4NeutronHPCaptureFS:
Public Member Functions | |
G4NeutronHPCaptureFS () | |
~G4NeutronHPCaptureFS () | |
void | Init (G4double A, G4double Z, G4int M, G4String &dirName, G4String &aFSType) |
G4HadFinalState * | ApplyYourself (const G4HadProjectile &theTrack) |
G4NeutronHPFinalState * | New () |
Definition at line 40 of file G4NeutronHPCaptureFS.hh.
G4NeutronHPCaptureFS::G4NeutronHPCaptureFS | ( | ) | [inline] |
Definition at line 44 of file G4NeutronHPCaptureFS.hh.
References G4NeutronHPFinalState::hasXsec.
00045 { 00046 hasXsec = false; 00047 hasExactMF6 = false; 00048 targetMass = 0; 00049 }
G4NeutronHPCaptureFS::~G4NeutronHPCaptureFS | ( | ) | [inline] |
G4HadFinalState * G4NeutronHPCaptureFS::ApplyYourself | ( | const G4HadProjectile & | theTrack | ) | [virtual] |
Reimplemented from G4NeutronHPFinalState.
Definition at line 47 of file G4NeutronHPCaptureFS.cc.
References G4HadFinalState::AddSecondary(), G4PhotonEvaporation::BreakItUp(), G4HadFinalState::Clear(), G4ParticleTable::FindIon(), G4UniformRand, G4Gamma::Gamma(), G4DynamicParticle::Get4Momentum(), G4HadProjectile::Get4Momentum(), G4Nucleus::GetBiasedThermalNucleus(), G4DynamicParticle::GetDefinition(), G4HadProjectile::GetDefinition(), G4ReactionProduct::GetKineticEnergy(), G4HadProjectile::GetKineticEnergy(), G4HadProjectile::GetMaterial(), G4ReactionProduct::GetMomentum(), G4NucleiProperties::GetNuclearMass(), G4HadFinalState::GetNumberOfSecondaries(), G4HadSecondary::GetParticle(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGMass(), G4NeutronHPPhotonDist::GetPhotons(), G4HadFinalState::GetSecondary(), G4Material::GetTemperature(), G4ReactionProduct::GetTotalEnergy(), G4NeutronHPFinalState::HasFSData(), G4ReactionProduct::Lorentz(), G4Neutron::Neutron(), G4INCL::Math::pi, G4NeutronHPEnAngCorrelation::Sample(), G4DynamicParticle::SetDefinition(), G4ReactionProduct::SetDefinition(), G4PhotonEvaporation::SetICM(), G4ReactionProduct::SetKineticEnergy(), G4DynamicParticle::SetMomentum(), G4ReactionProduct::SetMomentum(), G4NeutronHPEnAngCorrelation::SetNeutron(), G4HadFinalState::SetStatusChange(), G4NeutronHPEnAngCorrelation::SetTarget(), stopAndKill, G4NeutronHPFinalState::theBaseA, G4NeutronHPFinalState::theBaseZ, G4NeutronHPFinalState::theResult, and TRUE.
00048 { 00049 00050 G4int i; 00051 theResult.Clear(); 00052 // prepare neutron 00053 G4double eKinetic = theTrack.GetKineticEnergy(); 00054 const G4HadProjectile *incidentParticle = &theTrack; 00055 G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) ); 00056 theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() ); 00057 theNeutron.SetKineticEnergy( eKinetic ); 00058 00059 // prepare target 00060 G4ReactionProduct theTarget; 00061 G4Nucleus aNucleus; 00062 G4double eps = 0.0001; 00063 if(targetMass<500*MeV) 00064 targetMass = ( G4NucleiProperties::GetNuclearMass( static_cast<G4int>(theBaseA+eps) , static_cast<G4int>(theBaseZ+eps) )) / 00065 G4Neutron::Neutron()->GetPDGMass(); 00066 G4ThreeVector neutronVelocity = 1./G4Neutron::Neutron()->GetPDGMass()*theNeutron.GetMomentum(); 00067 G4double temperature = theTrack.GetMaterial()->GetTemperature(); 00068 theTarget = aNucleus.GetBiasedThermalNucleus(targetMass, neutronVelocity, temperature); 00069 00070 // go to nucleus rest system 00071 theNeutron.Lorentz(theNeutron, -1*theTarget); 00072 eKinetic = theNeutron.GetKineticEnergy(); 00073 00074 // dice the photons 00075 00076 G4ReactionProductVector * thePhotons = 0; 00077 if ( HasFSData() && !getenv ( "G4NEUTRONHP_USE_ONLY_PHOTONEVAPORATION" ) ) 00078 { 00079 //TK110430 00080 if ( hasExactMF6 ) 00081 { 00082 theMF6FinalState.SetTarget(theTarget); 00083 theMF6FinalState.SetNeutron(theNeutron); 00084 thePhotons = theMF6FinalState.Sample( eKinetic ); 00085 } 00086 else 00087 thePhotons = theFinalStatePhotons.GetPhotons(eKinetic); 00088 } 00089 else 00090 { 00091 G4ThreeVector aCMSMomentum = theNeutron.GetMomentum()+theTarget.GetMomentum(); 00092 G4LorentzVector p4(aCMSMomentum, theTarget.GetTotalEnergy() + theNeutron.GetTotalEnergy()); 00093 G4Fragment nucleus(static_cast<G4int>(theBaseA+1), static_cast<G4int>(theBaseZ) ,p4); 00094 G4PhotonEvaporation photonEvaporation; 00095 // T. K. add 00096 photonEvaporation.SetICM( TRUE ); 00097 G4FragmentVector* products = photonEvaporation.BreakItUp(nucleus); 00098 G4FragmentVector::iterator it; 00099 thePhotons = new G4ReactionProductVector; 00100 for(it=products->begin(); it!=products->end(); it++) 00101 { 00102 G4ReactionProduct * theOne = new G4ReactionProduct; 00103 // T. K. add 00104 if ( (*it)->GetParticleDefinition() != 0 ) 00105 theOne->SetDefinition( (*it)->GetParticleDefinition() ); 00106 else 00107 theOne->SetDefinition( G4Gamma::Gamma() ); // this definiion will be over writen 00108 00109 // T. K. comment out below line 00110 //theOne->SetDefinition( G4Gamma::Gamma() ); 00111 G4ParticleTable* theTable = G4ParticleTable::GetParticleTable(); 00112 if( (*it)->GetMomentum().mag() > 10*MeV ) theOne->SetDefinition( theTable->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ)) ); 00113 00114 //if ( (*i)->GetExcitationEnergy() > 0 ) 00115 if ( (*it)->GetExcitationEnergy() > 1.0e-2*eV ) 00116 { 00117 G4double ex = (*it)->GetExcitationEnergy(); 00118 G4ReactionProduct* aPhoton = new G4ReactionProduct; 00119 aPhoton->SetDefinition( G4Gamma::Gamma() ); 00120 aPhoton->SetMomentum( (*it)->GetMomentum().vect().unit() * ex ); 00121 //aPhoton->SetTotalEnergy( ex ); //will be calculated from momentum 00122 thePhotons->push_back(aPhoton); 00123 } 00124 00125 theOne->SetMomentum( (*it)->GetMomentum().vect() * ( (*it)->GetMomentum().t() - (*it)->GetExcitationEnergy() ) / (*it)->GetMomentum().t() ) ; 00126 //theOne->SetTotalEnergy( (*i)->GetMomentum().t() - (*i)->GetExcitationEnergy() ); //will be calculated from momentum 00127 thePhotons->push_back(theOne); 00128 delete *it; 00129 } 00130 delete products; 00131 } 00132 00133 00134 00135 // add them to the final state 00136 00137 G4int nPhotons = 0; 00138 if(thePhotons!=0) nPhotons=thePhotons->size(); 00139 00140 00141 //110527TKDB Unused codes, Detected by gcc4.6 compiler 00142 //G4int nParticles = nPhotons; 00143 //if(1==nPhotons) nParticles = 2; 00144 00145 //Make at least one photon 00146 //101203 TK 00147 if ( nPhotons == 0 ) 00148 { 00149 G4ReactionProduct * theOne = new G4ReactionProduct; 00150 theOne->SetDefinition( G4Gamma::Gamma() ); 00151 G4double theta = pi*G4UniformRand(); 00152 G4double phi = twopi*G4UniformRand(); 00153 G4double sinth = std::sin(theta); 00154 G4ThreeVector direction( sinth*std::cos(phi), sinth*std::sin(phi), std::cos(theta) ); 00155 theOne->SetMomentum( direction ) ; 00156 thePhotons->push_back(theOne); 00157 nPhotons++; // 0 -> 1 00158 } 00159 //One photon case: energy set to Q-value 00160 //101203 TK 00161 //if ( nPhotons == 1 ) 00162 if ( nPhotons == 1 && thePhotons->operator[](0)->GetDefinition()->GetBaryonNumber() == 0 ) 00163 { 00164 G4ThreeVector direction = thePhotons->operator[](0)->GetMomentum().unit(); 00165 G4double Q = G4ParticleTable::GetParticleTable()->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA), 0, static_cast<G4int>(theBaseZ))->GetPDGMass() + G4Neutron::Neutron()->GetPDGMass() 00166 - G4ParticleTable::GetParticleTable()->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ))->GetPDGMass(); 00167 thePhotons->operator[](0)->SetMomentum( Q*direction ); 00168 } 00169 // 00170 00171 // back to lab system 00172 for(i=0; i<nPhotons; i++) 00173 { 00174 thePhotons->operator[](i)->Lorentz(*(thePhotons->operator[](i)), theTarget); 00175 } 00176 00177 // Recoil, if only one gamma 00178 //if (1==nPhotons) 00179 if ( nPhotons == 1 && thePhotons->operator[](0)->GetDefinition()->GetBaryonNumber() == 0 ) 00180 { 00181 G4DynamicParticle * theOne = new G4DynamicParticle; 00182 G4ParticleDefinition * aRecoil = G4ParticleTable::GetParticleTable() 00183 ->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ)); 00184 theOne->SetDefinition(aRecoil); 00185 // Now energy; 00186 // Can be done slightly better @ 00187 G4ThreeVector aMomentum = theTrack.Get4Momentum().vect() 00188 +theTarget.GetMomentum() 00189 -thePhotons->operator[](0)->GetMomentum(); 00190 00191 G4ThreeVector theMomUnit = aMomentum.unit(); 00192 G4double aKinEnergy = theTrack.GetKineticEnergy() 00193 +theTarget.GetKineticEnergy(); // gammas come from Q-value 00194 G4double theResMass = aRecoil->GetPDGMass(); 00195 G4double theResE = aRecoil->GetPDGMass()+aKinEnergy; 00196 G4double theAbsMom = std::sqrt(theResE*theResE - theResMass*theResMass); 00197 G4ThreeVector theMomentum = theAbsMom*theMomUnit; 00198 theOne->SetMomentum(theMomentum); 00199 theResult.AddSecondary(theOne); 00200 } 00201 00202 // Now fill in the gammas. 00203 for(i=0; i<nPhotons; i++) 00204 { 00205 // back to lab system 00206 G4DynamicParticle * theOne = new G4DynamicParticle; 00207 theOne->SetDefinition(thePhotons->operator[](i)->GetDefinition()); 00208 theOne->SetMomentum(thePhotons->operator[](i)->GetMomentum()); 00209 theResult.AddSecondary(theOne); 00210 delete thePhotons->operator[](i); 00211 } 00212 delete thePhotons; 00213 00214 //101203TK 00215 G4bool residual = false; 00216 G4ParticleDefinition * aRecoil = G4ParticleTable::GetParticleTable() 00217 ->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ)); 00218 for ( G4int j = 0 ; j != theResult.GetNumberOfSecondaries() ; j++ ) 00219 { 00220 if ( theResult.GetSecondary(j)->GetParticle()->GetDefinition() == aRecoil ) residual = true; 00221 } 00222 00223 if ( residual == false ) 00224 { 00225 //G4ParticleDefinition * aRecoil = G4ParticleTable::GetParticleTable() 00226 // ->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ)); 00227 G4int nNonZero = 0; 00228 G4LorentzVector p_photons(0,0,0,0); 00229 for ( G4int j = 0 ; j != theResult.GetNumberOfSecondaries() ; j++ ) 00230 { 00231 p_photons += theResult.GetSecondary(j)->GetParticle()->Get4Momentum(); 00232 // To many 0 momentum photons -> Check PhotonDist 00233 if ( theResult.GetSecondary(j)->GetParticle()->Get4Momentum().e() > 0 ) nNonZero++; 00234 } 00235 00236 // Can we include kinetic energy here? 00237 G4double deltaE = ( theTrack.Get4Momentum().e() + theTarget.GetTotalEnergy() ) 00238 - ( p_photons.e() + aRecoil->GetPDGMass() ); 00239 00240 //Add photons 00241 if ( nPhotons - nNonZero > 0 ) 00242 { 00243 //G4cout << "TKDB G4NeutronHPCaptureFS::ApplyYourself we will create additional " << nPhotons - nNonZero << " photons" << G4endl; 00244 std::vector<G4double> vRand; 00245 vRand.push_back( 0.0 ); 00246 for ( G4int j = 0 ; j != nPhotons - nNonZero - 1 ; j++ ) 00247 { 00248 vRand.push_back( G4UniformRand() ); 00249 } 00250 vRand.push_back( 1.0 ); 00251 std::sort( vRand.begin(), vRand.end() ); 00252 00253 std::vector<G4double> vEPhoton; 00254 for ( G4int j = 0 ; j < (G4int)vRand.size() - 1 ; j++ ) 00255 { 00256 vEPhoton.push_back( deltaE * ( vRand[j+1] - vRand[j] ) ); 00257 } 00258 std::sort( vEPhoton.begin(), vEPhoton.end() ); 00259 00260 for ( G4int j = 0 ; j < nPhotons - nNonZero - 1 ; j++ ) 00261 { 00262 //Isotopic in LAB OK? 00263 G4double theta = pi*G4UniformRand(); 00264 G4double phi = twopi*G4UniformRand(); 00265 G4double sinth = std::sin(theta); 00266 G4double en = vEPhoton[j]; 00267 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) ); 00268 00269 p_photons += G4LorentzVector ( tempVector, tempVector.mag() ); 00270 G4DynamicParticle * theOne = new G4DynamicParticle; 00271 theOne->SetDefinition( G4Gamma::Gamma() ); 00272 theOne->SetMomentum( tempVector ); 00273 theResult.AddSecondary(theOne); 00274 } 00275 00276 // Add last photon 00277 G4DynamicParticle * theOne = new G4DynamicParticle; 00278 theOne->SetDefinition( G4Gamma::Gamma() ); 00279 // For better momentum conservation 00280 G4ThreeVector lastPhoton = -p_photons.vect().unit()*vEPhoton.back(); 00281 p_photons += G4LorentzVector( lastPhoton , lastPhoton.mag() ); 00282 theOne->SetMomentum( lastPhoton ); 00283 theResult.AddSecondary(theOne); 00284 } 00285 00286 //Add residual 00287 G4DynamicParticle * theOne = new G4DynamicParticle; 00288 G4ThreeVector aMomentum = theTrack.Get4Momentum().vect() + theTarget.GetMomentum() 00289 - p_photons.vect(); 00290 theOne->SetDefinition(aRecoil); 00291 theOne->SetMomentum( aMomentum ); 00292 theResult.AddSecondary(theOne); 00293 00294 } 00295 //101203TK END 00296 00297 // clean up the primary neutron 00298 theResult.SetStatusChange(stopAndKill); 00299 return &theResult; 00300 }
void G4NeutronHPCaptureFS::Init | ( | G4double | A, | |
G4double | Z, | |||
G4int | M, | |||
G4String & | dirName, | |||
G4String & | aFSType | |||
) | [virtual] |
Implements G4NeutronHPFinalState.
Definition at line 303 of file G4NeutronHPCaptureFS.cc.
References G4NeutronHPDataUsed::GetName(), G4NeutronHPNames::GetName(), G4NeutronHPPhotonDist::GetTargetMass(), G4NeutronHPFinalState::hasAnyData, G4NeutronHPFinalState::hasFSData, G4NeutronHPFinalState::hasXsec, G4NeutronHPEnAngCorrelation::Init(), G4NeutronHPPhotonDist::InitAngular(), G4NeutronHPPhotonDist::InitEnergies(), G4NeutronHPPhotonDist::InitMean(), G4NeutronHPFinalState::SetAZMs(), G4NeutronHPFinalState::theBaseA, and G4NeutronHPFinalState::theBaseZ.
00304 { 00305 00306 //TK110430 BEGIN 00307 std::stringstream ss; 00308 ss << static_cast<G4int>(Z); 00309 G4String sZ; 00310 ss >> sZ; 00311 ss.clear(); 00312 ss << static_cast<G4int>(A); 00313 G4String sA; 00314 ss >> sA; 00315 00316 ss.clear(); 00317 G4String sM; 00318 if ( M > 0 ) 00319 { 00320 ss << "m"; 00321 ss << M; 00322 ss >> sM; 00323 ss.clear(); 00324 } 00325 00326 G4String element_name = theNames.GetName( static_cast<G4int>(Z)-1 ); 00327 G4String filenameMF6 = dirName+"/FSMF6/"+sZ+"_"+sA+sM+"_"+element_name; 00328 std::ifstream dummyIFS(filenameMF6, std::ios::in); 00329 if ( dummyIFS.good() == true ) hasExactMF6=true; 00330 00331 //TK110430 Just for checking 00332 //ENDF-VII.0 no case (check done at 110430 00333 /* 00334 if ( hasExactMF6 == true ) 00335 { 00336 G4String filename = dirName+"FS/"+sZ+"_"+sA+"_"+element_name; 00337 std::ifstream dummyIFS(filename, std::ios::in); 00338 if ( dummyIFS.good() == true ) 00339 { 00340 G4cout << "TKDB Capture Both FS and FSMF6 are exist for Z = " << sZ << ", A = " << sA << G4endl;; 00341 } 00342 } 00343 */ 00344 00345 //TK110430 Only use MF6MT102 which has exactly same A and Z 00346 //Even _nat_ do not select and there is no _nat_ case in ENDF-VII.0 00347 if ( hasExactMF6 == true ) 00348 { 00349 std::ifstream theData(filenameMF6, std::ios::in); 00350 theMF6FinalState.Init(theData); 00351 theData.close(); 00352 return; 00353 } 00354 //TK110430 END 00355 00356 00357 G4String tString = "/FS"; 00358 G4bool dbool; 00359 G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, dirName, tString, dbool); 00360 00361 G4String filename = aFile.GetName(); 00362 SetAZMs( A, Z, M, aFile ); 00363 //theBaseA = A; 00364 //theBaseZ = G4int(Z+.5); 00365 if(!dbool || ( Z<2.5 && ( std::abs(theBaseZ - Z)>0.0001 || std::abs(theBaseA - A)>0.0001))) 00366 { 00367 hasAnyData = false; 00368 hasFSData = false; 00369 hasXsec = false; 00370 return; 00371 } 00372 std::ifstream theData(filename, std::ios::in); 00373 00374 hasFSData = theFinalStatePhotons.InitMean(theData); 00375 if(hasFSData) 00376 { 00377 targetMass = theFinalStatePhotons.GetTargetMass(); 00378 theFinalStatePhotons.InitAngular(theData); 00379 theFinalStatePhotons.InitEnergies(theData); 00380 } 00381 theData.close(); 00382 }
G4NeutronHPFinalState* G4NeutronHPCaptureFS::New | ( | ) | [inline, virtual] |
Implements G4NeutronHPFinalState.
Definition at line 57 of file G4NeutronHPCaptureFS.hh.
00058 { 00059 G4NeutronHPCaptureFS * theNew = new G4NeutronHPCaptureFS; 00060 return theNew; 00061 }