#include <G4NeutronHPPhotonDist.hh>
Public Member Functions | |
G4NeutronHPPhotonDist () | |
~G4NeutronHPPhotonDist () | |
G4bool | InitMean (std::ifstream &aDataFile) |
void | InitAngular (std::ifstream &aDataFile) |
void | InitEnergies (std::ifstream &aDataFile) |
void | InitPartials (std::ifstream &aDataFile) |
G4ReactionProductVector * | GetPhotons (G4double anEnergy) |
G4double | GetTargetMass () |
G4bool | NeedsCascade () |
G4double | GetLevelEnergy () |
Definition at line 55 of file G4NeutronHPPhotonDist.hh.
G4NeutronHPPhotonDist::G4NeutronHPPhotonDist | ( | ) | [inline] |
Definition at line 59 of file G4NeutronHPPhotonDist.hh.
00060 : repFlag( 0 ) 00061 , targetMass( 0.0 ) 00062 , nDiscrete( 0 ) 00063 , isoFlag( 0 ) 00064 , tabulationType( 0 ) 00065 , nDiscrete2( 0 ) 00066 , nIso( 0 ) 00067 , nPartials( 0 ) 00068 , theInternalConversionFlag( 0 ) 00069 , nGammaEnergies( 0 ) 00070 , theBaseEnergy( 0.0 ) 00071 { 00072 00073 disType = 0; 00074 energy = 0; 00075 theYield = 0; 00076 thePartialXsec = 0; 00077 isPrimary = 0; 00078 theShells = 0; 00079 theGammas = 0; 00080 nNeu = 0; 00081 theLegendre = 0; 00082 theAngular = 0; 00083 distribution = 0; 00084 probs = 0; 00085 partials = 0; 00086 actualMult = 0; 00087 00088 theLevelEnergies = 0; 00089 theTransitionProbabilities = 0; 00090 thePhotonTransitionFraction = 0; 00091 00092 }
G4NeutronHPPhotonDist::~G4NeutronHPPhotonDist | ( | ) | [inline] |
Definition at line 94 of file G4NeutronHPPhotonDist.hh.
00095 { 00096 delete [] disType; 00097 delete [] energy; 00098 delete [] theYield; 00099 delete [] thePartialXsec; 00100 delete [] isPrimary; 00101 delete [] theShells; 00102 delete [] theGammas; 00103 delete [] nNeu; 00104 delete [] theAngular; 00105 delete [] distribution; 00106 delete [] probs; 00107 00108 if ( theLegendre != NULL ) 00109 { 00110 for ( G4int i = 0 ; i < (nDiscrete2-nIso) ; i++ ) 00111 if ( theLegendre[i] != NULL ) delete[] theLegendre[i]; 00112 00113 delete [] theLegendre; 00114 } 00115 00116 if ( partials != 0 ) 00117 { 00118 for ( G4int i = 0 ; i < nPartials ; i++ ) 00119 { delete partials[i]; } 00120 00121 delete [] partials; 00122 } 00123 00124 delete [] actualMult; 00125 00126 // delete theLevelEnergies; 00127 // delete theTransitionProbabilities; 00128 // delete thePhotonTransitionFraction; 00129 // TKDB 00130 delete [] theLevelEnergies; 00131 delete [] theTransitionProbabilities; 00132 delete [] thePhotonTransitionFraction; 00133 }
G4double G4NeutronHPPhotonDist::GetLevelEnergy | ( | ) | [inline] |
Definition at line 149 of file G4NeutronHPPhotonDist.hh.
Referenced by G4NeutronHPInelasticCompFS::CompositeApply().
G4ReactionProductVector * G4NeutronHPPhotonDist::GetPhotons | ( | G4double | anEnergy | ) |
Definition at line 281 of file G4NeutronHPPhotonDist.cc.
References DBL_MAX, G4Electron::Electron(), G4cout, G4endl, G4Poisson(), G4UniformRand, G4Gamma::Gamma(), G4NeutronHPAngularP::GetCosTh(), G4ParticleDefinition::GetPDGMass(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), G4NeutronHPVector::GetVectorLength(), G4NeutronHPVector::GetX(), G4NeutronHPPartial::GetY(), G4INCL::Math::pi, G4NeutronHPVector::Sample(), G4NeutronHPLegendreStore::SampleMax(), G4NeutronHPLegendreStore::SetCoeff(), G4ReactionProduct::SetDefinition(), G4ReactionProduct::SetMomentum(), and G4ReactionProduct::SetTotalEnergy().
Referenced by G4NeutronHPCaptureFS::ApplyYourself(), G4NeutronHPInelasticBaseFS::BaseApply(), G4NeutronHPInelasticCompFS::CompositeApply(), and G4NeutronHPFSFissionFS::GetPhotons().
00282 { 00283 00284 //G4cout << "G4NeutronHPPhotonDist::GetPhotons repFlag " << repFlag << G4endl; 00285 // the partial cross-section case is not in this yet. @@@@ << 070601 TK add partial 00286 G4int i, ii, iii; 00287 G4int nSecondaries = 0; 00288 G4ReactionProductVector * thePhotons = new G4ReactionProductVector; 00289 if(repFlag==1) 00290 { 00291 G4double current=0; 00292 for(i=0; i<nDiscrete; i++) 00293 { 00294 current = theYield[i].GetY(anEnergy); 00295 actualMult[i] = G4Poisson(current); // max cut-off still missing @@@ 00296 if(nDiscrete==1&¤t<1.0001) 00297 { 00298 actualMult[i] = static_cast<G4int>(current); 00299 if(current<1) 00300 { 00301 actualMult[i] = 0; 00302 if(G4UniformRand()<current) actualMult[i] = 1; 00303 } 00304 } 00305 nSecondaries += actualMult[i]; 00306 } 00307 //G4cout << "nSecondaries " << nSecondaries << " anEnergy " << anEnergy/eV << G4endl; 00308 for(i=0;i<nSecondaries;i++) 00309 { 00310 G4ReactionProduct * theOne = new G4ReactionProduct; 00311 theOne->SetDefinition(G4Gamma::Gamma()); 00312 thePhotons->push_back(theOne); 00313 } 00314 G4int count=0; 00315 00316 /* 00317 G4double totalCascadeEnergy = 0.; 00318 G4double lastCascadeEnergy = 0.; 00319 G4double eGamm = 0; 00320 G4int maxEnergyIndex = 0; 00321 */ 00322 //Gcout << "nDiscrete " << nDiscrete << " nPartials " << nPartials << G4endl; 00323 //3456 00324 if ( nDiscrete == 1 && nPartials == 1 ) 00325 { 00326 if ( actualMult[ 0 ] > 0 ) 00327 { 00328 if ( disType[0] == 1 ) // continuum 00329 { 00330 00331 /* 00332 for(ii=0; ii< actualMult[0]; ii++) 00333 { 00334 00335 G4double sum=0, run=0; 00336 for(iii=0; iii<nPartials; iii++) sum+=probs[iii].GetY(anEnergy); 00337 G4double random = G4UniformRand(); 00338 G4int theP = 0; 00339 for(iii=0; iii<nPartials; iii++) 00340 { 00341 run+=probs[iii].GetY(anEnergy); 00342 theP = iii; 00343 if(random<run/sum) break; 00344 } 00345 if(theP==nPartials) theP=nPartials-1; // das sortiert J aus. 00346 sum=0; 00347 G4NeutronHPVector * temp; 00348 temp = partials[theP]->GetY(anEnergy); //@@@ look at, seems fishy 00349 // Looking for TotalCascdeEnergy or LastMaxEnergy 00350 if (ii == 0) 00351 { 00352 maxEnergyIndex = temp->GetVectorLength()-1; 00353 totalCascadeEnergy = temp->GetX(maxEnergyIndex); 00354 lastCascadeEnergy = totalCascadeEnergy; 00355 } 00356 lastCascadeEnergy -= eGamm; 00357 if (ii != actualMult[i]-1) eGamm = temp->SampleWithMax(lastCascadeEnergy); 00358 else eGamm = lastCascadeEnergy; 00359 thePhotons->operator[](count)->SetKineticEnergy(eGamm); 00360 delete temp; 00361 00362 } 00363 */ 00364 G4NeutronHPVector * temp; 00365 temp = partials[ 0 ]->GetY(anEnergy); //@@@ look at, seems fishy 00366 G4double maximumE = temp->GetX( temp->GetVectorLength()-1 ); // This is an assumption. 00367 00368 //G4cout << "start " << actualMult[ 0 ] << " maximumE " << maximumE/eV << G4endl; 00369 00370 std::vector< G4double > photons_e_best( actualMult[ 0 ] , 0.0 ); 00371 G4double best = DBL_MAX; 00372 G4int maxTry = 1000; 00373 for ( G4int j = 0 ; j < maxTry ; j++ ) 00374 { 00375 std::vector< G4double > photons_e( actualMult[ 0 ] , 0.0 ); 00376 for ( std::vector< G4double >::iterator 00377 it = photons_e.begin() ; it < photons_e.end() ; it++ ) 00378 { 00379 *it = temp->Sample(); 00380 } 00381 if ( std::accumulate( photons_e.begin() , photons_e.end() , 0.0 ) > maximumE ) 00382 { 00383 if ( std::accumulate( photons_e.begin() , photons_e.end() , 0.0 ) < best ) 00384 photons_e_best = photons_e; 00385 continue; 00386 } 00387 else 00388 { 00389 for ( std::vector< G4double >::iterator 00390 it = photons_e.begin() ; it < photons_e.end() ; it++ ) 00391 { 00392 thePhotons->operator[](count)->SetKineticEnergy( *it ); 00393 } 00394 //G4cout << "OK " << actualMult[0] << " j " << j << " total photons E " 00395 // << std::accumulate( photons_e.begin() , photons_e.end() , 0.0 )/eV << " ratio " << std::accumulate( photons_e.begin() , photons_e.end() , 0.0 ) / maximumE 00396 // << G4endl; 00397 00398 break; 00399 } 00400 G4cout << "NeutronHPPhotonDist could not find fitted energy set for multiplicity of " << actualMult[0] << "." << G4endl; 00401 G4cout << "NeutronHPPhotonDist will use the best set." << G4endl; 00402 for ( std::vector< G4double >::iterator 00403 it = photons_e_best.begin() ; it < photons_e_best.end() ; it++ ) 00404 { 00405 thePhotons->operator[](count)->SetKineticEnergy( *it ); 00406 } 00407 //G4cout << "Not Good " << actualMult[0] << " j " << j << " total photons E " 00408 // << best/eV << " ratio " << best / maximumE 00409 // << G4endl; 00410 } 00411 // TKDB 00412 delete temp; 00413 } 00414 else // discrete 00415 { 00416 thePhotons->operator[](count)->SetKineticEnergy(energy[i]); 00417 } 00418 count++; 00419 if(count > nSecondaries) throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist::GetPhotons inconsistancy"); 00420 } 00421 00422 } 00423 else 00424 { 00425 for(i=0; i<nDiscrete; i++) 00426 { 00427 for(ii=0; ii< actualMult[i]; ii++) 00428 { 00429 if(disType[i]==1) // continuum 00430 { 00431 G4double sum=0, run=0; 00432 for(iii=0; iii<nPartials; iii++) sum+=probs[iii].GetY(anEnergy); 00433 G4double random = G4UniformRand(); 00434 G4int theP = 0; 00435 for(iii=0; iii<nPartials; iii++) 00436 { 00437 run+=probs[iii].GetY(anEnergy); 00438 theP = iii; 00439 if(random<run/sum) break; 00440 } 00441 if(theP==nPartials) theP=nPartials-1; // das sortiert J aus. 00442 sum=0; 00443 G4NeutronHPVector * temp; 00444 temp = partials[theP]->GetY(anEnergy); //@@@ look at, seems fishy 00445 G4double eGamm = temp->Sample(); 00446 thePhotons->operator[](count)->SetKineticEnergy(eGamm); 00447 delete temp; 00448 } 00449 else // discrete 00450 { 00451 thePhotons->operator[](count)->SetKineticEnergy(energy[i]); 00452 } 00453 count++; 00454 if(count > nSecondaries) throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist::GetPhotons inconsistancy"); 00455 } 00456 } 00457 } 00458 // now do the angular distributions... 00459 if( isoFlag == 1) 00460 { 00461 for (i=0; i< nSecondaries; i++) 00462 { 00463 G4double costheta = 2.*G4UniformRand()-1; 00464 G4double theta = std::acos(costheta); 00465 G4double phi = twopi*G4UniformRand(); 00466 G4double sinth = std::sin(theta); 00467 G4double en = thePhotons->operator[](i)->GetTotalEnergy(); 00468 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) ); 00469 thePhotons->operator[](i)->SetMomentum( temp ) ; 00470 // G4cout << "Isotropic distribution in PhotonDist"<<temp<<G4endl; 00471 } 00472 } 00473 else 00474 { 00475 for(i=0; i<nSecondaries; i++) 00476 { 00477 G4double currentEnergy = thePhotons->operator[](i)->GetTotalEnergy(); 00478 for(ii=0; ii<nDiscrete2; ii++) 00479 { 00480 if (std::abs(currentEnergy-theGammas[ii])<0.1*keV) break; 00481 } 00482 if(ii==nDiscrete2) ii--; // fix for what seems an (file12 vs file 14) inconsistancy found in the ENDF 7N14 data. @@ 00483 if(ii<nIso) 00484 { 00485 // isotropic distribution 00486 G4double theta = pi*G4UniformRand(); 00487 G4double phi = twopi*G4UniformRand(); 00488 G4double sinth = std::sin(theta); 00489 G4double en = thePhotons->operator[](i)->GetTotalEnergy(); 00490 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) ); 00491 thePhotons->operator[](i)->SetMomentum( tempVector ) ; 00492 } 00493 else if(tabulationType==1) 00494 { 00495 // legendre polynomials 00496 G4int it(0); 00497 for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy 00498 { 00499 it = iii; 00500 if(theLegendre[ii-nIso][iii].GetEnergy()>anEnergy) 00501 break; 00502 } 00503 G4NeutronHPLegendreStore aStore(2); 00504 aStore.SetCoeff(1, &(theLegendre[ii-nIso][it])); 00505 //aStore.SetCoeff(0, &(theLegendre[ii-nIso][it-1])); 00506 //TKDB 110512 00507 if ( it > 0 ) 00508 { 00509 aStore.SetCoeff(0, &(theLegendre[ii-nIso][it-1])); 00510 } 00511 else 00512 { 00513 aStore.SetCoeff(0, &(theLegendre[ii-nIso][it])); 00514 } 00515 G4double cosTh = aStore.SampleMax(anEnergy); 00516 G4double theta = std::acos(cosTh); 00517 G4double phi = twopi*G4UniformRand(); 00518 G4double sinth = std::sin(theta); 00519 G4double en = thePhotons->operator[](i)->GetTotalEnergy(); 00520 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) ); 00521 thePhotons->operator[](i)->SetMomentum( tempVector ) ; 00522 } 00523 else 00524 { 00525 // tabulation of probabilities. 00526 G4int it(0); 00527 for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy 00528 { 00529 it = iii; 00530 if(theAngular[ii-nIso][iii].GetEnergy()>anEnergy) 00531 break; 00532 } 00533 G4double costh = theAngular[ii-nIso][it].GetCosTh(); // no interpolation yet @@ 00534 G4double theta = std::acos(costh); 00535 G4double phi = twopi*G4UniformRand(); 00536 G4double sinth = std::sin(theta); 00537 G4double en = thePhotons->operator[](i)->GetTotalEnergy(); 00538 G4ThreeVector tmpVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costh ); 00539 thePhotons->operator[](i)->SetMomentum( tmpVector ) ; 00540 } 00541 } 00542 } 00543 } 00544 else if(repFlag == 2) 00545 { 00546 G4double * running = new G4double[nGammaEnergies]; 00547 running[0]=theTransitionProbabilities[0]; 00548 //G4int i; //declaration at 284th 00549 for(i=1; i<nGammaEnergies; i++) 00550 { 00551 running[i]=running[i-1]+theTransitionProbabilities[i]; 00552 } 00553 G4double random = G4UniformRand(); 00554 G4int it=0; 00555 for(i=0; i<nGammaEnergies; i++) 00556 { 00557 it = i; 00558 if(random < running[i]/running[nGammaEnergies-1]) break; 00559 } 00560 delete [] running; 00561 G4double totalEnergy = theBaseEnergy - theLevelEnergies[it]; 00562 G4ReactionProduct * theOne = new G4ReactionProduct; 00563 theOne->SetDefinition(G4Gamma::Gamma()); 00564 random = G4UniformRand(); 00565 if(theInternalConversionFlag==2 && random>thePhotonTransitionFraction[it]) 00566 { 00567 theOne->SetDefinition(G4Electron::Electron()); 00568 //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009 00569 //But never enter at least with G4NDL3.13 00570 totalEnergy += G4Electron::Electron()->GetPDGMass(); //proposed correction: add this line for electron 00571 } 00572 theOne->SetTotalEnergy(totalEnergy); 00573 if( isoFlag == 1 ) 00574 { 00575 G4double costheta = 2.*G4UniformRand()-1; 00576 G4double theta = std::acos(costheta); 00577 G4double phi = twopi*G4UniformRand(); 00578 G4double sinth = std::sin(theta); 00579 //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009 00580 //G4double en = theOne->GetTotalEnergy(); 00581 G4double en = theOne->GetTotalMomentum(); 00582 //But never cause real effect at least with G4NDL3.13 TK 00583 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) ); 00584 theOne->SetMomentum( temp ) ; 00585 } 00586 else 00587 { 00588 G4double currentEnergy = theOne->GetTotalEnergy(); 00589 for(ii=0; ii<nDiscrete2; ii++) 00590 { 00591 if (std::abs(currentEnergy-theGammas[ii])<0.1*keV) break; 00592 } 00593 if(ii==nDiscrete2) ii--; // fix for what seems an (file12 vs file 14) inconsistancy found in the ENDF 7N14 data. @@ 00594 if(ii<nIso) 00595 { 00596 //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009 00597 // isotropic distribution 00598 //G4double theta = pi*G4UniformRand(); 00599 G4double theta = std::acos(2.*G4UniformRand()-1.); 00600 //But this is alos never cause real effect at least with G4NDL3.13 TK not repFlag == 2 AND isoFlag != 1 00601 G4double phi = twopi*G4UniformRand(); 00602 G4double sinth = std::sin(theta); 00603 //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009 00604 //G4double en = theOne->GetTotalEnergy(); 00605 G4double en = theOne->GetTotalMomentum(); 00606 //But never cause real effect at least with G4NDL3.13 TK 00607 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) ); 00608 theOne->SetMomentum( tempVector ) ; 00609 } 00610 else if(tabulationType==1) 00611 { 00612 // legendre polynomials 00613 G4int itt(0); 00614 for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy 00615 { 00616 itt = iii; 00617 if(theLegendre[ii-nIso][iii].GetEnergy()>anEnergy) 00618 break; 00619 } 00620 G4NeutronHPLegendreStore aStore(2); 00621 aStore.SetCoeff(1, &(theLegendre[ii-nIso][itt])); 00622 //aStore.SetCoeff(0, &(theLegendre[ii-nIso][it-1])); 00623 //TKDB 110512 00624 if ( itt > 0 ) 00625 { 00626 aStore.SetCoeff(0, &(theLegendre[ii-nIso][itt-1])); 00627 } 00628 else 00629 { 00630 aStore.SetCoeff(0, &(theLegendre[ii-nIso][itt])); 00631 } 00632 G4double cosTh = aStore.SampleMax(anEnergy); 00633 G4double theta = std::acos(cosTh); 00634 G4double phi = twopi*G4UniformRand(); 00635 G4double sinth = std::sin(theta); 00636 //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009 00637 //G4double en = theOne->GetTotalEnergy(); 00638 G4double en = theOne->GetTotalMomentum(); 00639 //But never cause real effect at least with G4NDL3.13 TK 00640 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) ); 00641 theOne->SetMomentum( tempVector ) ; 00642 } 00643 else 00644 { 00645 // tabulation of probabilities. 00646 G4int itt(0); 00647 for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy 00648 { 00649 itt = iii; 00650 if(theAngular[ii-nIso][iii].GetEnergy()>anEnergy) 00651 break; 00652 } 00653 G4double costh = theAngular[ii-nIso][itt].GetCosTh(); // no interpolation yet @@ 00654 G4double theta = std::acos(costh); 00655 G4double phi = twopi*G4UniformRand(); 00656 G4double sinth = std::sin(theta); 00657 //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009 00658 //G4double en = theOne->GetTotalEnergy(); 00659 G4double en = theOne->GetTotalMomentum(); 00660 //But never cause real effect at least with G4NDL3.13 TK 00661 G4ThreeVector tmpVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costh ); 00662 theOne->SetMomentum( tmpVector ) ; 00663 } 00664 } 00665 thePhotons->push_back(theOne); 00666 } 00667 else if( repFlag==0 ) 00668 { 00669 00670 // TK add 00671 if ( thePartialXsec == 0 ) 00672 { 00673 //G4cout << "repFlag is 0, but no PartialXsec data" << G4endl; 00674 //G4cout << "This is not support yet." << G4endl; 00675 return thePhotons; 00676 } 00677 00678 // Partial Case 00679 00680 G4ReactionProduct * theOne = new G4ReactionProduct; 00681 theOne->SetDefinition( G4Gamma::Gamma() ); 00682 thePhotons->push_back( theOne ); 00683 00684 // Energy 00685 00686 //G4cout << "Partial Case nDiscrete " << nDiscrete << G4endl; 00687 G4double sum = 0.0; 00688 std::vector < G4double > dif( nDiscrete , 0.0 ); 00689 for ( G4int j = 0 ; j < nDiscrete ; j++ ) 00690 { 00691 G4double x = thePartialXsec[ j ].GetXsec( anEnergy ); // x in barn 00692 if ( x > 0 ) 00693 { 00694 sum += x; 00695 } 00696 dif [ j ] = sum; 00697 //G4cout << "j " << j << ", x " << x << ", dif " << dif [ j ] << G4endl; 00698 } 00699 00700 G4double rand = G4UniformRand(); 00701 00702 G4int iphoton = 0; 00703 for ( G4int j = 0 ; j < nDiscrete ; j++ ) 00704 { 00705 G4double y = rand*sum; 00706 if ( dif [ j ] > y ) 00707 { 00708 iphoton = j; 00709 break; 00710 } 00711 } 00712 //G4cout << "iphoton " << iphoton << G4endl; 00713 //G4cout << "photon energy " << theGammas[ iphoton ] /eV << G4endl; 00714 00715 // Angle 00716 G4double cosTheta = 0.0; // mu 00717 00718 if ( isoFlag == 1 ) 00719 { 00720 00721 // Isotropic Case 00722 00723 cosTheta = 2.*G4UniformRand()-1; 00724 00725 } 00726 else 00727 { 00728 00729 if ( iphoton < nIso ) 00730 { 00731 00732 // still Isotropic 00733 00734 cosTheta = 2.*G4UniformRand()-1; 00735 00736 } 00737 else 00738 { 00739 00740 //G4cout << "Not Isotropic and isoFlag " << isoFlag << G4endl; 00741 //G4cout << "tabulationType " << tabulationType << G4endl; 00742 //G4cout << "nDiscrete2 " << nDiscrete2 << G4endl; 00743 //G4cout << "nIso " << nIso << G4endl; 00744 //G4cout << "size of nNeu " << nDiscrete2-nIso << G4endl; 00745 //G4cout << "nNeu[iphoton-nIso] " << nNeu[iphoton-nIso] << G4endl; 00746 00747 if ( tabulationType == 1 ) 00748 { 00749 // legendre polynomials 00750 00751 G4int iangle = 0; 00752 for ( G4int j = 0 ; j < nNeu [ iphoton - nIso ] ; j++ ) 00753 { 00754 iangle = j; 00755 if ( theLegendre[ iphoton - nIso ][ j ].GetEnergy() > anEnergy ) break; 00756 } 00757 00758 G4NeutronHPLegendreStore aStore( 2 ); 00759 aStore.SetCoeff( 1 , &( theLegendre[ iphoton - nIso ][ iangle ] ) ); 00760 aStore.SetCoeff( 0 , &( theLegendre[ iphoton - nIso ][ iangle - 1 ] ) ); 00761 00762 cosTheta = aStore.SampleMax( anEnergy ); 00763 00764 } 00765 else if ( tabulationType == 2 ) 00766 { 00767 00768 // tabulation of probabilities. 00769 00770 G4int iangle = 0; 00771 for ( G4int j = 0 ; j < nNeu [ iphoton - nIso ] ; j++ ) 00772 { 00773 iangle = j; 00774 if ( theAngular[ iphoton - nIso ][ j ].GetEnergy() > anEnergy ) break; 00775 } 00776 00777 cosTheta = theAngular[iphoton-nIso][ iangle ].GetCosTh(); // no interpolation yet @@ 00778 00779 } 00780 } 00781 } 00782 00783 // Set 00784 G4double phi = twopi*G4UniformRand(); 00785 G4double theta = std::acos( cosTheta ); 00786 G4double sinTheta = std::sin( theta ); 00787 00788 G4double photonE = theGammas[ iphoton ]; 00789 G4ThreeVector direction ( sinTheta*std::cos( phi ) , sinTheta * std::sin( phi ) , cosTheta ); 00790 G4ThreeVector photonP = photonE * direction; 00791 thePhotons->operator[]( 0 )->SetMomentum( photonP ) ; 00792 00793 } 00794 else 00795 { 00796 delete thePhotons; 00797 thePhotons = 0; // no gamma data available; some work needed @@@@@@@ 00798 } 00799 return thePhotons; 00800 }
G4double G4NeutronHPPhotonDist::GetTargetMass | ( | ) | [inline] |
Definition at line 145 of file G4NeutronHPPhotonDist.hh.
Referenced by G4NeutronHPCaptureFS::Init().
void G4NeutronHPPhotonDist::InitAngular | ( | std::ifstream & | aDataFile | ) |
Definition at line 121 of file G4NeutronHPPhotonDist.cc.
References G4cout, G4endl, and G4InterpolationManager::Init().
Referenced by G4NeutronHPInelasticCompFS::Init(), G4NeutronHPInelasticBaseFS::Init(), G4NeutronHPFSFissionFS::Init(), G4NeutronHPCaptureFS::Init(), and G4FissionLibrary::Init().
00122 { 00123 00124 G4int i, ii; 00125 //angular distributions 00126 aDataFile >> isoFlag; 00127 if (isoFlag != 1) 00128 { 00129 if ( repFlag == 2 ) G4cout << "G4NeutronHPPhotonDist: repFlag == 2 && isoFlag != 1 is unexpected! If you use G4ND3.x, then please report to Geant4 Hyper News. Thanks." << G4endl; 00130 aDataFile >> tabulationType >> nDiscrete2 >> nIso; 00131 //080731 00132 if ( theGammas != NULL && nDiscrete2 != nDiscrete ) G4cout << "080731c G4NeutronHPPhotonDist nDiscrete2 != nDiscrete, It looks like something wrong in your NDL files. Please update the latest. If you still have this messages after the update, then please report to Geant4 Hyper News." << G4endl; 00133 00134 // The order of cross section (InitPartials) and distribution (InitAngular here) data are different, we have to re-coordinate consistent data order. 00135 std::vector < G4double > vct_gammas_par; 00136 std::vector < G4double > vct_shells_par; 00137 std::vector < G4int > vct_primary_par; 00138 std::vector < G4int > vct_distype_par; 00139 std::vector < G4NeutronHPVector* > vct_pXS_par; 00140 if ( theGammas != NULL ) 00141 { 00142 //copy the cross section data 00143 for ( i = 0 ; i < nDiscrete ; i++ ) 00144 { 00145 vct_gammas_par.push_back( theGammas[ i ] ); 00146 vct_shells_par.push_back( theShells[ i ] ); 00147 vct_primary_par.push_back( isPrimary[ i ] ); 00148 vct_distype_par.push_back( disType[ i ] ); 00149 G4NeutronHPVector* hpv = new G4NeutronHPVector; 00150 *hpv = thePartialXsec[ i ]; 00151 vct_pXS_par.push_back( hpv ); 00152 } 00153 } 00154 if ( theGammas == NULL ) theGammas = new G4double[nDiscrete2]; 00155 if ( theShells == NULL ) theShells = new G4double[nDiscrete2]; 00156 //080731 00157 00158 for (i=0; i< nIso; i++) // isotropic photons 00159 { 00160 aDataFile >> theGammas[i] >> theShells[i]; 00161 theGammas[i]*=eV; 00162 theShells[i]*=eV; 00163 } 00164 nNeu = new G4int [nDiscrete2-nIso]; 00165 if(tabulationType==1)theLegendre=new G4NeutronHPLegendreTable *[nDiscrete2-nIso]; 00166 if(tabulationType==2)theAngular =new G4NeutronHPAngularP *[nDiscrete2-nIso]; 00167 for(i=nIso; i< nDiscrete2; i++) 00168 { 00169 if(tabulationType==1) 00170 { 00171 aDataFile >> theGammas[i] >> theShells[i] >> nNeu[i-nIso]; 00172 theGammas[i]*=eV; 00173 theShells[i]*=eV; 00174 theLegendre[i-nIso]=new G4NeutronHPLegendreTable[nNeu[i-nIso]]; 00175 theLegendreManager.Init(aDataFile); 00176 for (ii=0; ii<nNeu[i-nIso]; ii++) 00177 { 00178 theLegendre[i-nIso][ii].Init(aDataFile); 00179 } 00180 } 00181 else if(tabulationType==2) 00182 { 00183 aDataFile >> theGammas[i] >> theShells[i] >> nNeu[i-nIso]; 00184 theGammas[i]*=eV; 00185 theShells[i]*=eV; 00186 theAngular[i-nIso]=new G4NeutronHPAngularP[nNeu[i-nIso]]; 00187 for (ii=0; ii<nNeu[i-nIso]; ii++) 00188 { 00189 theAngular[i-nIso][ii].Init(aDataFile); 00190 } 00191 } 00192 else 00193 { 00194 G4cout << "tabulation type: tabulationType"<<G4endl; 00195 throw G4HadronicException(__FILE__, __LINE__, "cannot deal with this tabulation type for angular distributions."); 00196 } 00197 } 00198 //080731 00199 if ( vct_gammas_par.size() > 0 ) 00200 { 00201 //Reordering cross section data to corrsponding distribution data 00202 for ( i = 0 ; i < nDiscrete ; i++ ) 00203 { 00204 for ( G4int j = 0 ; j < nDiscrete ; j++ ) 00205 { 00206 // Checking gamma and shell to identification 00207 if ( theGammas[ i ] == vct_gammas_par [ j ] && theShells [ i ] == vct_shells_par[ j ] ) 00208 { 00209 isPrimary [ i ] = vct_primary_par [ j ]; 00210 disType [ i ] = vct_distype_par [ j ]; 00211 thePartialXsec[ i ] = ( *( vct_pXS_par[ j ] ) ); 00212 } 00213 } 00214 } 00215 //Garbage collection 00216 for ( std::vector < G4NeutronHPVector* >::iterator 00217 it = vct_pXS_par.begin() ; it != vct_pXS_par.end() ; it++ ) 00218 { 00219 delete *it; 00220 } 00221 } 00222 //080731 00223 } 00224 }
void G4NeutronHPPhotonDist::InitEnergies | ( | std::ifstream & | aDataFile | ) |
Definition at line 227 of file G4NeutronHPPhotonDist.cc.
References G4NeutronHPPartial::Init(), G4NeutronHPVector::Init(), and G4NeutronHPPartial::InitInterpolation().
Referenced by G4NeutronHPInelasticCompFS::Init(), G4NeutronHPInelasticBaseFS::Init(), G4NeutronHPFSFissionFS::Init(), G4NeutronHPCaptureFS::Init(), and G4FissionLibrary::Init().
00228 { 00229 G4int i, energyDistributionsNeeded = 0; 00230 for (i=0; i<nDiscrete; i++) 00231 { 00232 if( disType[i]==1) energyDistributionsNeeded =1; 00233 } 00234 if(!energyDistributionsNeeded) return; 00235 aDataFile >> nPartials; 00236 distribution = new G4int[nPartials]; 00237 probs = new G4NeutronHPVector[nPartials]; 00238 partials = new G4NeutronHPPartial * [nPartials]; 00239 G4int nen; 00240 G4int dummy; 00241 for (i=0; i<nPartials; i++) 00242 { 00243 aDataFile >> dummy; 00244 probs[i].Init(aDataFile, eV); 00245 aDataFile >> nen; 00246 partials[i] = new G4NeutronHPPartial(nen); 00247 partials[i]->InitInterpolation(aDataFile); 00248 partials[i]->Init(aDataFile); 00249 } 00250 }
G4bool G4NeutronHPPhotonDist::InitMean | ( | std::ifstream & | aDataFile | ) |
Definition at line 54 of file G4NeutronHPPhotonDist.cc.
References G4cout, G4endl, and G4NeutronHPVector::Init().
Referenced by G4NeutronHPInelasticCompFS::Init(), G4NeutronHPInelasticBaseFS::Init(), G4NeutronHPFSFissionFS::Init(), G4NeutronHPCaptureFS::Init(), and G4FissionLibrary::Init().
00055 { 00056 00057 G4bool result = true; 00058 if(aDataFile >> repFlag) 00059 { 00060 00061 aDataFile >> targetMass; 00062 if(repFlag==1) 00063 { 00064 // multiplicities 00065 aDataFile >> nDiscrete; 00066 disType = new G4int[nDiscrete]; 00067 energy = new G4double[nDiscrete]; 00068 actualMult = new G4int[nDiscrete]; 00069 theYield = new G4NeutronHPVector[nDiscrete]; 00070 for (G4int i=0; i<nDiscrete; i++) 00071 { 00072 aDataFile >> disType[i]>>energy[i]; 00073 energy[i]*=eV; 00074 theYield[i].Init(aDataFile, eV); 00075 } 00076 } 00077 else if(repFlag == 2) 00078 { 00079 aDataFile >> theInternalConversionFlag; 00080 aDataFile >> theBaseEnergy; 00081 theBaseEnergy*=eV; 00082 aDataFile >> theInternalConversionFlag; 00083 // theInternalConversionFlag == 1 No IC, theInternalConversionFlag == 2 with IC 00084 aDataFile >> nGammaEnergies; 00085 theLevelEnergies = new G4double[nGammaEnergies]; 00086 theTransitionProbabilities = new G4double[nGammaEnergies]; 00087 if(theInternalConversionFlag == 2) thePhotonTransitionFraction = new G4double[nGammaEnergies]; 00088 for(G4int ii=0; ii<nGammaEnergies; ii++) 00089 { 00090 if(theInternalConversionFlag == 1) 00091 { 00092 aDataFile >> theLevelEnergies[ii] >> theTransitionProbabilities[ii]; 00093 theLevelEnergies[ii]*=eV; 00094 } 00095 else if(theInternalConversionFlag == 2) 00096 { 00097 aDataFile >> theLevelEnergies[ii] >> theTransitionProbabilities[ii] >> thePhotonTransitionFraction[ii]; 00098 theLevelEnergies[ii]*=eV; 00099 } 00100 else 00101 { 00102 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist: Unknown conversion flag"); 00103 } 00104 } 00105 // Note, that this is equivalent to using the 'Gamma' classes. 00106 // throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist: Transition probability array not sampled for the moment."); 00107 } 00108 else 00109 { 00110 G4cout << "Data representation in G4NeutronHPPhotonDist: "<<repFlag<<G4endl; 00111 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist: This data representation is not implemented."); 00112 } 00113 } 00114 else 00115 { 00116 result = false; 00117 } 00118 return result; 00119 }
void G4NeutronHPPhotonDist::InitPartials | ( | std::ifstream & | aDataFile | ) |
Definition at line 252 of file G4NeutronHPPhotonDist.cc.
References G4NeutronHPVector::Init().
Referenced by G4NeutronHPInelasticCompFS::Init(), and G4NeutronHPInelasticBaseFS::Init().
00253 { 00254 00255 //G4cout << "G4NeutronHPPhotonDist::InitPartials " << G4endl; 00256 aDataFile >> nDiscrete >> targetMass; 00257 if(nDiscrete != 1) 00258 { 00259 theTotalXsec.Init(aDataFile, eV); 00260 } 00261 G4int i; 00262 theGammas = new G4double[nDiscrete]; 00263 theShells = new G4double[nDiscrete]; 00264 isPrimary = new G4int[nDiscrete]; 00265 disType = new G4int[nDiscrete]; 00266 thePartialXsec = new G4NeutronHPVector[nDiscrete]; 00267 for(i=0; i<nDiscrete; i++) 00268 { 00269 aDataFile>>theGammas[i]>>theShells[i]>>isPrimary[i]>>disType[i]; 00270 theGammas[i]*=eV; 00271 theShells[i]*=eV; 00272 thePartialXsec[i].Init(aDataFile, eV); 00273 } 00274 00275 //G4cout << "G4NeutronHPPhotonDist::InitPartials Test " << G4endl; 00276 //G4cout << "G4NeutronHPPhotonDist::InitPartials nDiscrete " << nDiscrete << G4endl; 00277 //G4NeutronHPVector* aHP = new G4NeutronHPVector; 00278 //aHP->Check(1); 00279 }
G4bool G4NeutronHPPhotonDist::NeedsCascade | ( | ) | [inline] |