G4NeutronHPCaptureFS Class Reference

#include <G4NeutronHPCaptureFS.hh>

Inheritance diagram for G4NeutronHPCaptureFS:

G4NeutronHPFinalState

Public Member Functions

 G4NeutronHPCaptureFS ()
 ~G4NeutronHPCaptureFS ()
void Init (G4double A, G4double Z, G4int M, G4String &dirName, G4String &aFSType)
G4HadFinalStateApplyYourself (const G4HadProjectile &theTrack)
G4NeutronHPFinalStateNew ()

Detailed Description

Definition at line 40 of file G4NeutronHPCaptureFS.hh.


Constructor & Destructor Documentation

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]

Definition at line 51 of file G4NeutronHPCaptureFS.hh.

00052   {
00053   }


Member Function Documentation

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   }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:36 2013 for Geant4 by  doxygen 1.4.7