#include <G4NeutronHPInelasticCompFS.hh>
Inheritance diagram for G4NeutronHPInelasticCompFS:
Definition at line 42 of file G4NeutronHPInelasticCompFS.hh.
G4NeutronHPInelasticCompFS::G4NeutronHPInelasticCompFS | ( | ) | [inline] |
Definition at line 46 of file G4NeutronHPInelasticCompFS.hh.
References G4NeutronHPFinalState::hasXsec, LR, QI, theAngularDistribution, theEnergyAngData, theEnergyDistribution, theFinalStatePhotons, and theXsection.
00047 { 00048 00049 QI.resize(51); 00050 LR.resize(51); 00051 for(G4int i=0; i<51; i++) 00052 { 00053 hasXsec = true; 00054 theXsection[i] = 0; 00055 theEnergyDistribution[i] = 0; 00056 theAngularDistribution[i] = 0; 00057 theEnergyAngData[i] = 0; 00058 theFinalStatePhotons[i] = 0; 00059 QI[i]=0.0; 00060 LR[i]=0; 00061 } 00062 00063 }
virtual G4NeutronHPInelasticCompFS::~G4NeutronHPInelasticCompFS | ( | ) | [inline, virtual] |
Definition at line 64 of file G4NeutronHPInelasticCompFS.hh.
References theAngularDistribution, theEnergyAngData, theEnergyDistribution, theFinalStatePhotons, and theXsection.
00065 { 00066 for(G4int i=0; i<51; i++) 00067 { 00068 if(theXsection[i] != 0) delete theXsection[i]; 00069 if(theEnergyDistribution[i] != 0) delete theEnergyDistribution[i]; 00070 if(theAngularDistribution[i] != 0) delete theAngularDistribution[i]; 00071 if(theEnergyAngData[i] != 0) delete theEnergyAngData[i]; 00072 if(theFinalStatePhotons[i] != 0) delete theFinalStatePhotons[i]; 00073 } 00074 }
virtual G4HadFinalState* G4NeutronHPInelasticCompFS::ApplyYourself | ( | const G4HadProjectile & | theTrack | ) | [pure virtual] |
Reimplemented from G4NeutronHPFinalState.
Implemented in G4NeutronHPAInelasticFS, G4NeutronHPDInelasticFS, G4NeutronHPHe3InelasticFS, G4NeutronHPNInelasticFS, G4NeutronHPPInelasticFS, and G4NeutronHPTInelasticFS.
void G4NeutronHPInelasticCompFS::CompositeApply | ( | const G4HadProjectile & | theTrack, | |
G4ParticleDefinition * | aHadron | |||
) |
Definition at line 216 of file G4NeutronHPInelasticCompFS.cc.
References G4HadFinalState::AddSecondary(), G4NeutronHPFinalState::adjust_final_state(), G4HadFinalState::Clear(), G4UniformRand, G4Gamma::Gamma(), G4DynamicParticle::Get4Momentum(), G4HadProjectile::Get4Momentum(), G4ParticleDefinition::GetBaryonNumber(), G4Nucleus::GetBiasedThermalNucleus(), G4NeutronHPDeExGammas::GetDecayGammas(), G4ReactionProduct::GetDefinition(), G4HadProjectile::GetDefinition(), G4ParticleTable::GetIon(), G4ReactionProduct::GetKineticEnergy(), G4HadProjectile::GetKineticEnergy(), G4NeutronHPDeExGammas::GetLevel(), G4NeutronHPPhotonDist::GetLevelEnergy(), G4NeutronHPLevel::GetLevelEnergy(), G4NeutronHPDeExGammas::GetLevelEnergy(), G4ReactionProduct::GetMass(), G4HadProjectile::GetMaterial(), G4ReactionProduct::GetMomentum(), G4NucleiProperties::GetNuclearMass(), G4NeutronHPDeExGammas::GetNumberOfLevels(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4NeutronHPPhotonDist::GetPhotons(), G4Material::GetTemperature(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), InitDistributionInitialState(), G4ReactionProduct::Lorentz(), G4Neutron::Neutron(), QI, G4NeutronHPEnAngCorrelation::Sample(), G4NeutronHPEnergyDistribution::Sample(), G4NeutronHPAngular::SampleAndUpdate(), SelectExitChannel(), G4DynamicParticle::SetDefinition(), G4ReactionProduct::SetDefinition(), G4ReactionProduct::SetKineticEnergy(), G4DynamicParticle::SetMomentum(), G4ReactionProduct::SetMomentum(), G4HadFinalState::SetStatusChange(), stopAndKill, theAngularDistribution, G4NeutronHPFinalState::theBaseA, G4NeutronHPFinalState::theBaseZ, theEnergyAngData, theEnergyDistribution, theFinalStatePhotons, theGammas, G4NeutronHPFinalState::theResult, and theXsection.
Referenced by G4NeutronHPTInelasticFS::ApplyYourself(), G4NeutronHPPInelasticFS::ApplyYourself(), G4NeutronHPNInelasticFS::ApplyYourself(), G4NeutronHPHe3InelasticFS::ApplyYourself(), G4NeutronHPDInelasticFS::ApplyYourself(), and G4NeutronHPAInelasticFS::ApplyYourself().
00217 { 00218 00219 // prepare neutron 00220 theResult.Clear(); 00221 G4double eKinetic = theTrack.GetKineticEnergy(); 00222 const G4HadProjectile *incidentParticle = &theTrack; 00223 G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) ); 00224 theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() ); 00225 theNeutron.SetKineticEnergy( eKinetic ); 00226 00227 // prepare target 00228 G4int i; 00229 for(i=0; i<50; i++) 00230 { if(theXsection[i] != 0) { break; } } 00231 00232 G4double targetMass=0; 00233 G4double eps = 0.0001; 00234 targetMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps))) / 00235 G4Neutron::Neutron()->GetPDGMass(); 00236 // if(theEnergyAngData[i]!=0) 00237 // targetMass = theEnergyAngData[i]->GetTargetMass(); 00238 // else if(theAngularDistribution[i]!=0) 00239 // targetMass = theAngularDistribution[i]->GetTargetMass(); 00240 // else if(theFinalStatePhotons[50]!=0) 00241 // targetMass = theFinalStatePhotons[50]->GetTargetMass(); 00242 G4Nucleus aNucleus; 00243 G4ReactionProduct theTarget; 00244 G4ThreeVector neuVelo = (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum(); 00245 theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature()); 00246 00247 // prepare the residual mass 00248 G4double residualMass=0; 00249 G4double residualZ = theBaseZ - aDefinition->GetPDGCharge(); 00250 G4double residualA = theBaseA - aDefinition->GetBaryonNumber()+1; 00251 residualMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(residualA+eps), static_cast<G4int>(residualZ+eps)) ) / 00252 G4Neutron::Neutron()->GetPDGMass(); 00253 00254 // prepare energy in target rest frame 00255 G4ReactionProduct boosted; 00256 boosted.Lorentz(theNeutron, theTarget); 00257 eKinetic = boosted.GetKineticEnergy(); 00258 // G4double momentumInCMS = boosted.GetTotalMomentum(); 00259 00260 // select exit channel for composite FS class. 00261 G4int it = SelectExitChannel( eKinetic ); 00262 00263 // set target and neutron in the relevant exit channel 00264 InitDistributionInitialState(theNeutron, theTarget, it); 00265 00266 G4ReactionProductVector * thePhotons = 0; 00267 G4ReactionProductVector * theParticles = 0; 00268 G4ReactionProduct aHadron; 00269 aHadron.SetDefinition(aDefinition); // what if only cross-sections exist ==> Na 23 11 @@@@ 00270 G4double availableEnergy = theNeutron.GetKineticEnergy() + theNeutron.GetMass() - aHadron.GetMass() + 00271 (targetMass - residualMass)*G4Neutron::Neutron()->GetPDGMass(); 00272 //080730c 00273 if ( availableEnergy < 0 ) 00274 { 00275 //G4cout << "080730c Adjust availavleEnergy " << G4endl; 00276 availableEnergy = 0; 00277 } 00278 G4int nothingWasKnownOnHadron = 0; 00279 G4int dummy; 00280 G4double eGamm = 0; 00281 G4int iLevel=it-1; 00282 00283 // TK without photon has it = 0 00284 if( 50 == it ) 00285 { 00286 00287 // TK Excitation level is not determined 00288 iLevel=-1; 00289 aHadron.SetKineticEnergy(availableEnergy*residualMass*G4Neutron::Neutron()->GetPDGMass()/ 00290 (aHadron.GetMass()+residualMass*G4Neutron::Neutron()->GetPDGMass())); 00291 00292 //aHadron.SetMomentum(theNeutron.GetMomentum()*(1./theNeutron.GetTotalMomentum())* 00293 // std::sqrt(aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()- 00294 // aHadron.GetMass()*aHadron.GetMass())); 00295 00296 //TK add safty 100909 00297 G4double p2 = ( aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy() - aHadron.GetMass()*aHadron.GetMass() ); 00298 G4double p = 0.0; 00299 if ( p2 > 0.0 ) p = std::sqrt( p ); 00300 00301 aHadron.SetMomentum(theNeutron.GetMomentum()*(1./theNeutron.GetTotalMomentum())*p ); 00302 00303 } 00304 else 00305 { 00306 while( iLevel!=-1 && theGammas.GetLevel(iLevel) == 0 ) { iLevel--; } 00307 } 00308 00309 00310 if ( theAngularDistribution[it] != 0 ) // MF4 00311 { 00312 if(theEnergyDistribution[it]!=0) // MF5 00313 { 00314 //************************************************************ 00315 /* 00316 aHadron.SetKineticEnergy(theEnergyDistribution[it]->Sample(eKinetic, dummy)); 00317 G4double eSecN = aHadron.GetKineticEnergy(); 00318 */ 00319 //************************************************************ 00320 //EMendoza --> maximum allowable energy should be taken into account. 00321 G4double dqi = 0.0; 00322 if ( QI[it] < 0 || 849 < QI[it] ) dqi = QI[it]; //For backword compatibility QI introduced since G4NDL3.15 00323 G4double MaxEne=eKinetic+dqi; 00324 G4double eSecN; 00325 do{ 00326 eSecN=theEnergyDistribution[it]->Sample(eKinetic, dummy); 00327 }while(eSecN>MaxEne); 00328 aHadron.SetKineticEnergy(eSecN); 00329 //************************************************************ 00330 eGamm = eKinetic-eSecN; 00331 for(iLevel=theGammas.GetNumberOfLevels()-1; iLevel>=0; iLevel--) 00332 { 00333 if(theGammas.GetLevelEnergy(iLevel)<eGamm) break; 00334 } 00335 G4double random = 2*G4UniformRand(); 00336 iLevel+=G4int(random); 00337 if(iLevel>theGammas.GetNumberOfLevels()-1)iLevel = theGammas.GetNumberOfLevels()-1; 00338 } 00339 else 00340 { 00341 G4double eExcitation = 0; 00342 if(iLevel>=0) eExcitation = theGammas.GetLevel(iLevel)->GetLevelEnergy(); 00343 while (eKinetic-eExcitation < 0 && iLevel>0) 00344 { 00345 iLevel--; 00346 eExcitation = theGammas.GetLevel(iLevel)->GetLevelEnergy(); 00347 } 00348 //110610TK BEGIN 00349 //Use QI value for calculating excitation energy of residual. 00350 G4bool useQI=false; 00351 G4double dqi = QI[it]; 00352 if ( dqi < 0 || 849 < dqi ) useQI = true; //Former libraies does not have values of this range 00353 00354 if ( useQI ) 00355 { 00356 // QI introudced since G4NDL3.15 00357 eExcitation = -QI[it]; 00358 //Re-evluate iLevel based on this eExcitation 00359 iLevel = 0; 00360 G4bool find = false; 00361 G4int imaxEx = 0; 00362 while( !theGammas.GetLevel(iLevel+1) == 0 ) 00363 { 00364 G4double maxEx = 0.0; 00365 if ( maxEx < theGammas.GetLevel(iLevel)->GetLevelEnergy() ) 00366 { 00367 maxEx = theGammas.GetLevel(iLevel)->GetLevelEnergy(); 00368 imaxEx = iLevel; 00369 } 00370 if ( eExcitation < theGammas.GetLevel(iLevel)->GetLevelEnergy() ) 00371 { 00372 find = true; 00373 iLevel--; 00374 // very small eExcitation, iLevel becomes -1, this is protected below. 00375 if ( iLevel == -1 ) iLevel = 0; // But cause energy trouble. 00376 break; 00377 } 00378 iLevel++; 00379 } 00380 // In case, cannot find proper level, then use the maximum level. 00381 if ( !find ) iLevel = imaxEx; 00382 } 00383 //110610TK END 00384 00385 if(getenv("InelasticCompFSLogging") && eKinetic-eExcitation < 0) 00386 { 00387 throw G4HadronicException(__FILE__, __LINE__, "SEVERE: InelasticCompFS: Consistency of data not good enough, please file report"); 00388 } 00389 if(eKinetic-eExcitation < 0) eExcitation = 0; 00390 if(iLevel!= -1) aHadron.SetKineticEnergy(eKinetic - eExcitation); 00391 00392 } 00393 theAngularDistribution[it]->SampleAndUpdate(aHadron); 00394 00395 if( theFinalStatePhotons[it] == 0 ) 00396 { 00397 //G4cout << "110610 USE Gamma Level" << G4endl; 00398 // TK comment Most n,n* eneter to this 00399 thePhotons = theGammas.GetDecayGammas(iLevel); 00400 eGamm -= theGammas.GetLevelEnergy(iLevel); 00401 if(eGamm>0) // @ ok for now, but really needs an efficient way of correllated sampling @ 00402 { 00403 G4ReactionProduct * theRestEnergy = new G4ReactionProduct; 00404 theRestEnergy->SetDefinition(G4Gamma::Gamma()); 00405 theRestEnergy->SetKineticEnergy(eGamm); 00406 G4double costh = 2.*G4UniformRand()-1.; 00407 G4double phi = twopi*G4UniformRand(); 00408 theRestEnergy->SetMomentum(eGamm*std::sin(std::acos(costh))*std::cos(phi), 00409 eGamm*std::sin(std::acos(costh))*std::sin(phi), 00410 eGamm*costh); 00411 if(thePhotons == 0) { thePhotons = new G4ReactionProductVector; } 00412 thePhotons->push_back(theRestEnergy); 00413 } 00414 } 00415 } 00416 else if(theEnergyAngData[it] != 0) // MF6 00417 { 00418 theParticles = theEnergyAngData[it]->Sample(eKinetic); 00419 } 00420 else 00421 { 00422 // @@@ what to do, if we have photon data, but no info on the hadron itself 00423 nothingWasKnownOnHadron = 1; 00424 } 00425 00426 //G4cout << "theFinalStatePhotons it " << it << G4endl; 00427 //G4cout << "theFinalStatePhotons[it] " << theFinalStatePhotons[it] << G4endl; 00428 //G4cout << "theFinalStatePhotons it " << it << G4endl; 00429 //G4cout << "theFinalStatePhotons[it] " << theFinalStatePhotons[it] << G4endl; 00430 //G4cout << "thePhotons " << thePhotons << G4endl; 00431 00432 if ( theFinalStatePhotons[it] != 0 ) 00433 { 00434 // the photon distributions are in the Nucleus rest frame. 00435 // TK residual rest frame 00436 G4ReactionProduct boosted_tmp; 00437 boosted_tmp.Lorentz(theNeutron, theTarget); 00438 G4double anEnergy = boosted_tmp.GetKineticEnergy(); 00439 thePhotons = theFinalStatePhotons[it]->GetPhotons(anEnergy); 00440 G4double aBaseEnergy = theFinalStatePhotons[it]->GetLevelEnergy(); 00441 G4double testEnergy = 0; 00442 if(thePhotons!=0 && thePhotons->size()!=0) 00443 { aBaseEnergy-=thePhotons->operator[](0)->GetTotalEnergy(); } 00444 if(theFinalStatePhotons[it]->NeedsCascade()) 00445 { 00446 while(aBaseEnergy>0.01*keV) 00447 { 00448 // cascade down the levels 00449 G4bool foundMatchingLevel = false; 00450 G4int closest = 2; 00451 G4double deltaEold = -1; 00452 for(G4int j=1; j<it; j++) 00453 { 00454 if(theFinalStatePhotons[j]!=0) 00455 { 00456 testEnergy = theFinalStatePhotons[j]->GetLevelEnergy(); 00457 } 00458 else 00459 { 00460 testEnergy = 0; 00461 } 00462 G4double deltaE = std::abs(testEnergy-aBaseEnergy); 00463 if(deltaE<0.1*keV) 00464 { 00465 G4ReactionProductVector * theNext = 00466 theFinalStatePhotons[j]->GetPhotons(anEnergy); 00467 thePhotons->push_back(theNext->operator[](0)); 00468 aBaseEnergy = testEnergy-theNext->operator[](0)->GetTotalEnergy(); 00469 delete theNext; 00470 foundMatchingLevel = true; 00471 break; // ===> 00472 } 00473 if(theFinalStatePhotons[j]!=0 && ( deltaE<deltaEold||deltaEold<0.) ) 00474 { 00475 closest = j; 00476 deltaEold = deltaE; 00477 } 00478 } // <=== the break goes here. 00479 if(!foundMatchingLevel) 00480 { 00481 G4ReactionProductVector * theNext = 00482 theFinalStatePhotons[closest]->GetPhotons(anEnergy); 00483 thePhotons->push_back(theNext->operator[](0)); 00484 aBaseEnergy = aBaseEnergy-theNext->operator[](0)->GetTotalEnergy(); 00485 delete theNext; 00486 } 00487 } 00488 } 00489 } 00490 unsigned int i0; 00491 if(thePhotons!=0) 00492 { 00493 for(i0=0; i0<thePhotons->size(); i0++) 00494 { 00495 // back to lab 00496 thePhotons->operator[](i0)->Lorentz(*(thePhotons->operator[](i0)), -1.*theTarget); 00497 } 00498 } 00499 //G4cout << "nothingWasKnownOnHadron " << nothingWasKnownOnHadron << G4endl; 00500 if(nothingWasKnownOnHadron) 00501 { 00502 // TKDB 100405 00503 // In this case, hadron should be isotropic in CM 00504 // mu and p should be correlated 00505 // 00506 G4double totalPhotonEnergy = 0.0; 00507 if ( thePhotons != 0 ) 00508 { 00509 unsigned int nPhotons = thePhotons->size(); 00510 unsigned int ii0; 00511 for ( ii0=0; ii0<nPhotons; ii0++) 00512 { 00513 //thePhotons has energies at LAB system 00514 totalPhotonEnergy += thePhotons->operator[](ii0)->GetTotalEnergy(); 00515 } 00516 } 00517 00518 //isotropic distribution in CM 00519 G4double mu = 1.0 - 2 * G4UniformRand(); 00520 00521 // need momentums in target rest frame; 00522 G4LorentzVector target_in_LAB ( theTarget.GetMomentum() , theTarget.GetTotalEnergy() ); 00523 G4ThreeVector boostToTargetRest = -target_in_LAB.boostVector(); 00524 G4LorentzVector proj_in_LAB = incidentParticle->Get4Momentum(); 00525 00526 G4DynamicParticle* proj = new G4DynamicParticle( G4Neutron::Neutron() , proj_in_LAB.boost( boostToTargetRest ) ); 00527 G4DynamicParticle* targ = new G4DynamicParticle( G4ParticleTable::GetParticleTable()->GetIon ( (G4int)theBaseZ , (G4int)theBaseA , totalPhotonEnergy ) , G4ThreeVector(0) ); 00528 G4DynamicParticle* hadron = new G4DynamicParticle( aHadron.GetDefinition() , G4ThreeVector(0) ); // will be fill momentum 00529 00530 two_body_reaction ( proj , targ , hadron , mu ); 00531 00532 G4LorentzVector hadron_in_trag_rest = hadron->Get4Momentum(); 00533 G4LorentzVector hadron_in_LAB = hadron_in_trag_rest.boost ( -boostToTargetRest ); 00534 aHadron.SetMomentum( hadron_in_LAB.v() ); 00535 aHadron.SetKineticEnergy ( hadron_in_LAB.e() - hadron_in_LAB.m() ); 00536 00537 delete proj; 00538 delete targ; 00539 delete hadron; 00540 00541 //TKDB 100405 00542 /* 00543 G4double totalPhotonEnergy = 0; 00544 if(thePhotons!=0) 00545 { 00546 unsigned int nPhotons = thePhotons->size(); 00547 unsigned int i0; 00548 for(i0=0; i0<nPhotons; i0++) 00549 { 00550 totalPhotonEnergy += thePhotons->operator[](i0)->GetTotalEnergy(); 00551 } 00552 } 00553 availableEnergy -= totalPhotonEnergy; 00554 residualMass += totalPhotonEnergy/G4Neutron::Neutron()->GetPDGMass(); 00555 aHadron.SetKineticEnergy(availableEnergy*residualMass*G4Neutron::Neutron()->GetPDGMass()/ 00556 (aHadron.GetMass()+residualMass*G4Neutron::Neutron()->GetPDGMass())); 00557 G4double CosTheta = 1.0 - 2.0*G4UniformRand(); 00558 G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta); 00559 G4double Phi = twopi*G4UniformRand(); 00560 G4ThreeVector Vector(std::cos(Phi)*SinTheta, std::sin(Phi)*SinTheta, CosTheta); 00561 //aHadron.SetMomentum(Vector* std::sqrt(aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()- 00562 // aHadron.GetMass()*aHadron.GetMass())); 00563 G4double p2 = aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()- aHadron.GetMass()*aHadron.GetMass(); 00564 00565 G4double p = 0.0; 00566 if ( p2 > 0.0 ) 00567 p = std::sqrt ( p2 ); 00568 00569 aHadron.SetMomentum( Vector*p ); 00570 */ 00571 00572 } 00573 00574 // fill the result 00575 // Beware - the recoil is not necessarily in the particles... 00576 // Can be calculated from momentum conservation? 00577 // The idea is that the particles ar emitted forst, and the gammas only once the 00578 // recoil is on the residual; assumption is that gammas do not contribute to 00579 // the recoil. 00580 // This needs more design @@@ 00581 00582 G4int nSecondaries = 2; // the hadron and the recoil 00583 G4bool needsSeparateRecoil = false; 00584 G4int totalBaryonNumber = 0; 00585 G4int totalCharge = 0; 00586 G4ThreeVector totalMomentum(0); 00587 if(theParticles != 0) 00588 { 00589 nSecondaries = theParticles->size(); 00590 G4ParticleDefinition * aDef; 00591 unsigned int ii0; 00592 for(ii0=0; ii0<theParticles->size(); ii0++) 00593 { 00594 aDef = theParticles->operator[](ii0)->GetDefinition(); 00595 totalBaryonNumber+=aDef->GetBaryonNumber(); 00596 totalCharge+=G4int(aDef->GetPDGCharge()+eps); 00597 totalMomentum += theParticles->operator[](ii0)->GetMomentum(); 00598 } 00599 if(totalBaryonNumber!=G4int(theBaseA+eps+incidentParticle->GetDefinition()->GetBaryonNumber())) 00600 { 00601 needsSeparateRecoil = true; 00602 nSecondaries++; 00603 residualA = G4int(theBaseA+eps+incidentParticle->GetDefinition()->GetBaryonNumber() 00604 -totalBaryonNumber); 00605 residualZ = G4int(theBaseZ+eps+incidentParticle->GetDefinition()->GetPDGCharge() 00606 -totalCharge); 00607 } 00608 } 00609 00610 G4int nPhotons = 0; 00611 if(thePhotons!=0) { nPhotons = thePhotons->size(); } 00612 nSecondaries += nPhotons; 00613 00614 G4DynamicParticle * theSec; 00615 00616 if( theParticles==0 ) 00617 { 00618 theSec = new G4DynamicParticle; 00619 theSec->SetDefinition(aHadron.GetDefinition()); 00620 theSec->SetMomentum(aHadron.GetMomentum()); 00621 theResult.AddSecondary(theSec); 00622 00623 aHadron.Lorentz(aHadron, theTarget); 00624 G4ReactionProduct theResidual; 00625 theResidual.SetDefinition(G4ParticleTable::GetParticleTable() 00626 ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0)); 00627 theResidual.SetKineticEnergy(aHadron.GetKineticEnergy()*aHadron.GetMass()/theResidual.GetMass()); 00628 00629 //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #6 00630 //theResidual.SetMomentum(-1.*aHadron.GetMomentum()); 00631 G4ThreeVector incidentNeutronMomentum = theNeutron.GetMomentum(); 00632 theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum()); 00633 00634 theResidual.Lorentz(theResidual, -1.*theTarget); 00635 G4ThreeVector totalPhotonMomentum(0,0,0); 00636 if(thePhotons!=0) 00637 { 00638 for(i=0; i<nPhotons; i++) 00639 { 00640 totalPhotonMomentum += thePhotons->operator[](i)->GetMomentum(); 00641 } 00642 } 00643 theSec = new G4DynamicParticle; 00644 theSec->SetDefinition(theResidual.GetDefinition()); 00645 theSec->SetMomentum(theResidual.GetMomentum()-totalPhotonMomentum); 00646 theResult.AddSecondary(theSec); 00647 } 00648 else 00649 { 00650 for(i0=0; i0<theParticles->size(); i0++) 00651 { 00652 theSec = new G4DynamicParticle; 00653 theSec->SetDefinition(theParticles->operator[](i0)->GetDefinition()); 00654 theSec->SetMomentum(theParticles->operator[](i0)->GetMomentum()); 00655 theResult.AddSecondary(theSec); 00656 delete theParticles->operator[](i0); 00657 } 00658 delete theParticles; 00659 if(needsSeparateRecoil && residualZ!=0) 00660 { 00661 G4ReactionProduct theResidual; 00662 theResidual.SetDefinition(G4ParticleTable::GetParticleTable() 00663 ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0)); 00664 G4double resiualKineticEnergy = theResidual.GetMass()*theResidual.GetMass(); 00665 resiualKineticEnergy += totalMomentum*totalMomentum; 00666 resiualKineticEnergy = std::sqrt(resiualKineticEnergy) - theResidual.GetMass(); 00667 // cout << "Kinetic energy of the residual = "<<resiualKineticEnergy<<endl; 00668 theResidual.SetKineticEnergy(resiualKineticEnergy); 00669 00670 //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4 00671 //theResidual.SetMomentum(-1.*totalMomentum); 00672 //G4ThreeVector incidentNeutronMomentum = theNeutron.GetMomentum(); 00673 //theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum()); 00674 //080717 TK Comment still do NOT include photon's mometum which produce by thePhotons 00675 theResidual.SetMomentum( theNeutron.GetMomentum() + theTarget.GetMomentum() - totalMomentum ); 00676 00677 theSec = new G4DynamicParticle; 00678 theSec->SetDefinition(theResidual.GetDefinition()); 00679 theSec->SetMomentum(theResidual.GetMomentum()); 00680 theResult.AddSecondary(theSec); 00681 } 00682 } 00683 if(thePhotons!=0) 00684 { 00685 for(i=0; i<nPhotons; i++) 00686 { 00687 theSec = new G4DynamicParticle; 00688 //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009 00689 //theSec->SetDefinition(G4Gamma::Gamma()); 00690 theSec->SetDefinition( thePhotons->operator[](i)->GetDefinition() ); 00691 //But never cause real effect at least with G4NDL3.13 TK 00692 theSec->SetMomentum(thePhotons->operator[](i)->GetMomentum()); 00693 theResult.AddSecondary(theSec); 00694 delete thePhotons->operator[](i); 00695 } 00696 // some garbage collection 00697 delete thePhotons; 00698 } 00699 00700 //080721 00701 G4ParticleDefinition* targ_pd = G4ParticleTable::GetParticleTable()->GetIon ( (G4int)theBaseZ , (G4int)theBaseA , 0.0 ); 00702 G4LorentzVector targ_4p_lab ( theTarget.GetMomentum() , std::sqrt( targ_pd->GetPDGMass()*targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2() ) ); 00703 G4LorentzVector proj_4p_lab = theTrack.Get4Momentum(); 00704 G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab; 00705 adjust_final_state ( init_4p_lab ); 00706 00707 // clean up the primary neutron 00708 theResult.SetStatusChange(stopAndKill); 00709 }
virtual G4NeutronHPVector* G4NeutronHPInelasticCompFS::GetXsec | ( | ) | [inline, virtual] |
Reimplemented from G4NeutronHPFinalState.
Definition at line 83 of file G4NeutronHPInelasticCompFS.hh.
References theXsection.
Referenced by SelectExitChannel().
00083 { return theXsection[50]; }
Reimplemented from G4NeutronHPFinalState.
Definition at line 79 of file G4NeutronHPInelasticCompFS.hh.
References theXsection.
00080 { 00081 return std::max(0., theXsection[50]->GetY(anEnergy)); 00082 }
void G4NeutronHPInelasticCompFS::Init | ( | G4double | A, | |
G4double | Z, | |||
G4int | M, | |||
G4String & | dirName, | |||
G4String & | aSFType | |||
) | [virtual] |
Implements G4NeutronHPFinalState.
Reimplemented in G4NeutronHPAInelasticFS, G4NeutronHPDInelasticFS, G4NeutronHPHe3InelasticFS, G4NeutronHPNInelasticFS, G4NeutronHPPInelasticFS, and G4NeutronHPTInelasticFS.
Definition at line 78 of file G4NeutronHPInelasticCompFS.cc.
References G4cout, G4endl, gammaPath, G4NeutronHPDataUsed::GetName(), G4NeutronHPNames::GetName(), G4NeutronHPFinalState::hasAnyData, G4NeutronHPFinalState::hasFSData, G4NeutronHPFinalState::hasXsec, G4NeutronHPEnAngCorrelation::Init(), G4NeutronHPEnergyDistribution::Init(), G4NeutronHPAngular::Init(), G4NeutronHPVector::Init(), G4NeutronHPPhotonDist::InitAngular(), G4NeutronHPPhotonDist::InitEnergies(), G4NeutronHPPhotonDist::InitMean(), G4NeutronHPPhotonDist::InitPartials(), LR, QI, G4NeutronHPFinalState::SetAZMs(), theAngularDistribution, theEnergyAngData, theEnergyDistribution, theFinalStatePhotons, G4NeutronHPFinalState::theNames, G4NeutronHPFinalState::theNDLDataA, G4NeutronHPFinalState::theNDLDataZ, and theXsection.
Referenced by G4NeutronHPTInelasticFS::Init(), G4NeutronHPPInelasticFS::Init(), G4NeutronHPNInelasticFS::Init(), G4NeutronHPHe3InelasticFS::Init(), G4NeutronHPDInelasticFS::Init(), and G4NeutronHPAInelasticFS::Init().
00079 { 00080 gammaPath = "/Inelastic/Gammas/"; 00081 if(!getenv("G4NEUTRONHPDATA")) 00082 throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files."); 00083 G4String tBase = getenv("G4NEUTRONHPDATA"); 00084 gammaPath = tBase+gammaPath; 00085 G4String tString = dirName; 00086 G4bool dbool; 00087 G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, tString, aFSType, dbool); 00088 G4String filename = aFile.GetName(); 00089 SetAZMs( A, Z, M, aFile ); 00090 //theBaseA = aFile.GetA(); 00091 //theBaseZ = aFile.GetZ(); 00092 //theNDLDataA = (int)aFile.GetA(); 00093 //theNDLDataZ = aFile.GetZ(); 00094 //if(!dbool || ( Z<2.5 && ( std::abs(theBaseZ - Z)>0.0001 || std::abs(theBaseA - A)>0.0001))) 00095 if ( !dbool || ( Z<2.5 && ( std::abs(theNDLDataZ - Z)>0.0001 || std::abs(theNDLDataA - A)>0.0001)) ) 00096 { 00097 if(getenv("NeutronHPNamesLogging")) G4cout << "Skipped = "<< filename <<" "<<A<<" "<<Z<<G4endl; 00098 hasAnyData = false; 00099 hasFSData = false; 00100 hasXsec = false; 00101 return; 00102 } 00103 //theBaseA = A; 00104 //theBaseZ = G4int(Z+.5); 00105 std::ifstream theData(filename, std::ios::in); 00106 if(!theData) 00107 { 00108 hasAnyData = false; 00109 hasFSData = false; 00110 hasXsec = false; 00111 theData.close(); 00112 return; 00113 } 00114 // here we go 00115 G4int infoType, dataType, dummy; 00116 G4int sfType, it; 00117 hasFSData = false; 00118 while (theData >> infoType) 00119 { 00120 hasFSData = true; 00121 theData >> dataType; 00122 theData >> sfType >> dummy; 00123 it = 50; 00124 if(sfType>=600||(sfType<100&&sfType>=50)) it = sfType%50; 00125 if(dataType==3) 00126 { 00127 //theData >> dummy >> dummy; 00128 //TK110430 00129 // QI and LR introudced since G4NDL3.15 00130 G4double dqi; 00131 G4int ilr; 00132 theData >> dqi >> ilr; 00133 00134 QI[ it ] = dqi*eV; 00135 LR[ it ] = ilr; 00136 theXsection[it] = new G4NeutronHPVector; 00137 G4int total; 00138 theData >> total; 00139 theXsection[it]->Init(theData, total, eV); 00140 //std::cout << theXsection[it]->GetXsec(1*MeV) << std::endl; 00141 } 00142 else if(dataType==4) 00143 { 00144 theAngularDistribution[it] = new G4NeutronHPAngular; 00145 theAngularDistribution[it]->Init(theData); 00146 } 00147 else if(dataType==5) 00148 { 00149 theEnergyDistribution[it] = new G4NeutronHPEnergyDistribution; 00150 theEnergyDistribution[it]->Init(theData); 00151 } 00152 else if(dataType==6) 00153 { 00154 theEnergyAngData[it] = new G4NeutronHPEnAngCorrelation; 00155 theEnergyAngData[it]->Init(theData); 00156 } 00157 else if(dataType==12) 00158 { 00159 theFinalStatePhotons[it] = new G4NeutronHPPhotonDist; 00160 theFinalStatePhotons[it]->InitMean(theData); 00161 } 00162 else if(dataType==13) 00163 { 00164 theFinalStatePhotons[it] = new G4NeutronHPPhotonDist; 00165 theFinalStatePhotons[it]->InitPartials(theData); 00166 } 00167 else if(dataType==14) 00168 { 00169 theFinalStatePhotons[it]->InitAngular(theData); 00170 } 00171 else if(dataType==15) 00172 { 00173 theFinalStatePhotons[it]->InitEnergies(theData); 00174 } 00175 else 00176 { 00177 throw G4HadronicException(__FILE__, __LINE__, "Data-type unknown to G4NeutronHPInelasticCompFS"); 00178 } 00179 } 00180 theData.close(); 00181 }
void G4NeutronHPInelasticCompFS::InitDistributionInitialState | ( | G4ReactionProduct & | aNeutron, | |
G4ReactionProduct & | aTarget, | |||
G4int | it | |||
) | [inline] |
Definition at line 86 of file G4NeutronHPInelasticCompFS.hh.
References G4NeutronHPEnAngCorrelation::SetNeutron(), G4NeutronHPAngular::SetNeutron(), G4NeutronHPEnAngCorrelation::SetTarget(), G4NeutronHPAngular::SetTarget(), theAngularDistribution, and theEnergyAngData.
Referenced by CompositeApply().
00089 { 00090 if(theAngularDistribution[it]!=0) 00091 { 00092 theAngularDistribution[it]->SetTarget(aTarget); 00093 theAngularDistribution[it]->SetNeutron(aNeutron); 00094 } 00095 if(theEnergyAngData[it]!=0) 00096 { 00097 theEnergyAngData[it]->SetTarget(aTarget); 00098 theEnergyAngData[it]->SetNeutron(aNeutron); 00099 } 00100 }
Definition at line 57 of file G4NeutronHPInelasticCompFS.cc.
References gammaPath, G4NeutronHPDeExGammas::Init(), and theGammas.
Referenced by G4NeutronHPTInelasticFS::Init(), G4NeutronHPPInelasticFS::Init(), G4NeutronHPNInelasticFS::Init(), G4NeutronHPHe3InelasticFS::Init(), G4NeutronHPDInelasticFS::Init(), and G4NeutronHPAInelasticFS::Init().
00058 { 00059 // char the[100] = {""}; 00060 // std::ostrstream ost(the, 100, std::ios::out); 00061 // ost <<gammaPath<<"z"<<ZR<<".a"<<AR; 00062 // G4String * aName = new G4String(the); 00063 // std::ifstream from(*aName, std::ios::in); 00064 00065 std::ostringstream ost; 00066 ost <<gammaPath<<"z"<<ZR<<".a"<<AR; 00067 G4String aName = ost.str(); 00068 std::ifstream from(aName, std::ios::in); 00069 00070 if(!from) return; // no data found for this isotope 00071 // std::ifstream theGammaData(*aName, std::ios::in); 00072 std::ifstream theGammaData(aName, std::ios::in); 00073 00074 theGammas.Init(theGammaData); 00075 // delete aName; 00076 }
virtual G4NeutronHPFinalState* G4NeutronHPInelasticCompFS::New | ( | ) | [pure virtual] |
Implements G4NeutronHPFinalState.
Implemented in G4NeutronHPAInelasticFS, G4NeutronHPDInelasticFS, G4NeutronHPHe3InelasticFS, G4NeutronHPNInelasticFS, G4NeutronHPPInelasticFS, and G4NeutronHPTInelasticFS.
Definition at line 183 of file G4NeutronHPInelasticCompFS.cc.
References G4UniformRand, GetXsec(), and theXsection.
Referenced by CompositeApply().
00184 { 00185 00186 // it = 0 has without Photon 00187 G4double running[50]; 00188 running[0] = 0; 00189 unsigned int i; 00190 for(i=0; i<50; i++) 00191 { 00192 if(i!=0) running[i]=running[i-1]; 00193 if(theXsection[i] != 0) 00194 { 00195 running[i] += std::max(0., theXsection[i]->GetXsec(eKinetic)); 00196 } 00197 } 00198 G4double random = G4UniformRand(); 00199 G4double sum = running[49]; 00200 G4int it = 50; 00201 if(0!=sum) 00202 { 00203 G4int i0; 00204 for(i0=0; i0<50; i0++) 00205 { 00206 it = i0; 00207 if(random < running[i0]/sum) break; 00208 } 00209 } 00210 //debug: it = 1; 00211 return it; 00212 }
G4String G4NeutronHPInelasticCompFS::gammaPath [protected] |
Definition at line 112 of file G4NeutronHPInelasticCompFS.hh.
Referenced by Init(), and InitGammas().
std::vector<G4int > G4NeutronHPInelasticCompFS::LR [protected] |
Definition at line 119 of file G4NeutronHPInelasticCompFS.hh.
Referenced by G4NeutronHPInelasticCompFS(), and Init().
std::vector< G4double > G4NeutronHPInelasticCompFS::QI [protected] |
Definition at line 118 of file G4NeutronHPInelasticCompFS.hh.
Referenced by CompositeApply(), G4NeutronHPInelasticCompFS(), and Init().
G4NeutronHPAngular* G4NeutronHPInelasticCompFS::theAngularDistribution[51] [protected] |
Definition at line 106 of file G4NeutronHPInelasticCompFS.hh.
Referenced by CompositeApply(), G4NeutronHPInelasticCompFS(), Init(), InitDistributionInitialState(), and ~G4NeutronHPInelasticCompFS().
G4double G4NeutronHPInelasticCompFS::theCurrentA [protected] |
Definition at line 114 of file G4NeutronHPInelasticCompFS.hh.
G4double G4NeutronHPInelasticCompFS::theCurrentZ [protected] |
Definition at line 115 of file G4NeutronHPInelasticCompFS.hh.
G4NeutronHPEnAngCorrelation* G4NeutronHPInelasticCompFS::theEnergyAngData[51] [protected] |
Definition at line 107 of file G4NeutronHPInelasticCompFS.hh.
Referenced by CompositeApply(), G4NeutronHPInelasticCompFS(), Init(), InitDistributionInitialState(), and ~G4NeutronHPInelasticCompFS().
Definition at line 105 of file G4NeutronHPInelasticCompFS.hh.
Referenced by CompositeApply(), G4NeutronHPInelasticCompFS(), Init(), and ~G4NeutronHPInelasticCompFS().
G4NeutronHPPhotonDist* G4NeutronHPInelasticCompFS::theFinalStatePhotons[51] [protected] |
Definition at line 109 of file G4NeutronHPInelasticCompFS.hh.
Referenced by CompositeApply(), G4NeutronHPInelasticCompFS(), Init(), and ~G4NeutronHPInelasticCompFS().
Definition at line 111 of file G4NeutronHPInelasticCompFS.hh.
Referenced by CompositeApply(), and InitGammas().
G4NeutronHPVector* G4NeutronHPInelasticCompFS::theXsection[51] [protected] |
Definition at line 104 of file G4NeutronHPInelasticCompFS.hh.
Referenced by CompositeApply(), G4NeutronHPInelasticCompFS(), GetXsec(), Init(), SelectExitChannel(), and ~G4NeutronHPInelasticCompFS().