G4LENDModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // Class Description
00027 // Final state production model for a LEND (Low Energy Nuclear Data) 
00028 // LEND is Geant4 interface for GIDI (General Interaction Data Interface) 
00029 // which gives a discription of nuclear and atomic reactions, such as
00030 //    Binary collision cross sections
00031 //    Particle number multiplicity distributions of reaction products
00032 //    Energy and angular distributions of reaction products
00033 //    Derived calculational constants
00034 // GIDI is developped at Lawrence Livermore National Laboratory
00035 // Class Description - End
00036 
00037 // 071025 First implementation done by T. Koi (SLAC/SCCS)
00038 // 101118 Name modifications for release T. Koi (SLAC/PPA)
00039 
00040 #include "G4LENDModel.hh"
00041 #include "G4PhysicalConstants.hh"
00042 #include "G4SystemOfUnits.hh"
00043 #include "G4NistManager.hh"
00044 
00045 G4LENDModel::G4LENDModel( G4String name )
00046 :G4HadronicInteraction( name )
00047 {
00048 
00049    SetMinEnergy( 0.*eV );
00050    SetMaxEnergy( 20.*MeV );
00051 
00052    default_evaluation = "endl99"; 
00053    allow_nat = false;
00054    allow_any = false;
00055 
00056    lend_manager = G4LENDManager::GetInstance();  
00057 
00058 }
00059 
00060 G4LENDModel::~G4LENDModel()
00061 {
00062    for ( std::map< G4int , G4LENDUsedTarget* >::iterator 
00063          it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
00064    { 
00065       delete it->second;  
00066    }
00067 }
00068 
00069 
00070 void G4LENDModel::recreate_used_target_map()
00071 {
00072 
00073    for ( std::map< G4int , G4LENDUsedTarget* >::iterator 
00074          it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
00075    { 
00076       delete it->second;  
00077    }
00078    usedTarget_map.clear();
00079 
00080    create_used_target_map();
00081 
00082 }
00083 
00084 
00085 
00086 void G4LENDModel::create_used_target_map()
00087 {
00088 
00089    lend_manager->RequestChangeOfVerboseLevel( verboseLevel );
00090 
00091    size_t numberOfElements = G4Element::GetNumberOfElements();
00092    static const G4ElementTable* theElementTable = G4Element::GetElementTable();
00093 
00094    for ( size_t i = 0 ; i < numberOfElements ; ++i )
00095    {
00096 
00097       const G4Element* anElement = (*theElementTable)[i];
00098       G4int numberOfIsotope = anElement->GetNumberOfIsotopes(); 
00099 
00100       if ( numberOfIsotope > 0 )
00101       {
00102       // User Defined Abundances   
00103          for ( G4int i_iso = 0 ; i_iso < numberOfIsotope ; i_iso++ )
00104          {
00105             G4int iZ = anElement->GetIsotope( i_iso )->GetZ();
00106             G4int iA = anElement->GetIsotope( i_iso )->GetN();
00107 
00108             G4LENDUsedTarget* aTarget = new G4LENDUsedTarget ( proj , default_evaluation , iZ , iA );  
00109             if ( allow_nat == true ) aTarget->AllowNat();
00110             if ( allow_any == true ) aTarget->AllowAny();
00111             usedTarget_map.insert( std::pair< G4int , G4LENDUsedTarget* > ( lend_manager->GetNucleusEncoding( iZ , iA ) , aTarget ) );
00112          }
00113       }
00114       else
00115       {
00116       // Natural Abundances   
00117          G4NistElementBuilder* nistElementBuild = lend_manager->GetNistElementBuilder();
00118          G4int iZ = int ( anElement->GetZ() );
00119          //G4cout << nistElementBuild->GetNumberOfNistIsotopes( int ( anElement->GetZ() ) ) << G4endl;
00120          G4int numberOfNistIso = nistElementBuild->GetNumberOfNistIsotopes( int ( anElement->GetZ() ) ); 
00121 
00122          for ( G4int ii = 0 ; ii < numberOfNistIso ; ii++ )
00123          {
00124             //G4cout << nistElementBuild->GetIsotopeAbundance( iZ , nistElementBuild->GetNistFirstIsotopeN( iZ ) + i ) << G4endl;
00125             if ( nistElementBuild->GetIsotopeAbundance( iZ , nistElementBuild->GetNistFirstIsotopeN( iZ ) + ii ) > 0 )
00126             {
00127                G4int iMass = nistElementBuild->GetNistFirstIsotopeN( iZ ) + ii;  
00128                //G4cout << iZ << " " << nistElementBuild->GetNistFirstIsotopeN( iZ ) + i << " " << nistElementBuild->GetIsotopeAbundance ( iZ , iMass ) << G4endl;  
00129 
00130                G4LENDUsedTarget* aTarget = new G4LENDUsedTarget ( proj , default_evaluation , iZ , iMass );  
00131                if ( allow_nat == true ) aTarget->AllowNat();
00132                if ( allow_any == true ) aTarget->AllowAny();
00133                usedTarget_map.insert( std::pair< G4int , G4LENDUsedTarget* > ( lend_manager->GetNucleusEncoding( iZ , iMass ) , aTarget ) );
00134 
00135             }
00136 
00137          }
00138 
00139       }
00140    }
00141 
00142 
00143 
00144    G4cout << "Dump UsedTarget for " << GetModelName() << G4endl;
00145    G4cout << "Requested Evaluation, Z , A -> Actual Evaluation, Z , A(0=Nat) , Pointer of Target" << G4endl;
00146    for ( std::map< G4int , G4LENDUsedTarget* >::iterator 
00147          it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
00148    {
00149       G4cout 
00150          << " " << it->second->GetWantedEvaluation() 
00151          << ", " << it->second->GetWantedZ() 
00152          << ", " << it->second->GetWantedA() 
00153          << " -> " << it->second->GetActualEvaluation() 
00154          << ", " << it->second->GetActualZ() 
00155          << ", " << it->second->GetActualA() 
00156          << ", " << it->second->GetTarget() 
00157          << G4endl; 
00158    } 
00159 
00160 }
00161   
00162 
00163   
00164 #include "G4ParticleTable.hh"
00165   
00166 G4HadFinalState * G4LENDModel::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& aTarg )
00167 {
00168 
00169    G4double temp = aTrack.GetMaterial()->GetTemperature();
00170 
00171    //G4int iZ = int ( aTarg.GetZ() );
00172    //G4int iA = int ( aTarg.GetN() );
00173    //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this) 
00174    G4int iZ = aTarg.GetZ_asInt();
00175    G4int iA = aTarg.GetA_asInt();
00176 
00177    G4double ke = aTrack.GetKineticEnergy();
00178 
00179    G4HadFinalState* theResult = new G4HadFinalState();
00180 
00181    G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA ) )->second->GetTarget();
00182 
00183    G4double aMu = aTarget->getElasticFinalState( ke*MeV, temp, NULL, NULL );
00184 
00185    G4double phi = twopi*G4UniformRand();
00186    G4double theta = std::acos( aMu );
00187    //G4double sinth = std::sin( theta );
00188 
00189    G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>( aTrack.GetDefinition() ) );
00190    theNeutron.SetMomentum( aTrack.Get4Momentum().vect() );
00191    theNeutron.SetKineticEnergy( ke );
00192 
00193    G4ReactionProduct theTarget( G4ParticleTable::GetParticleTable()->FindIon( iZ , iA , 0 , iZ ) );
00194 
00195    G4double mass = G4ParticleTable::GetParticleTable()->FindIon( iZ , iA , 0 , iZ )->GetPDGMass();
00196 
00197 // add Thermal motion 
00198    G4double kT = k_Boltzmann*temp;
00199    G4ThreeVector v ( G4RandGauss::shoot() * std::sqrt( kT*mass ) 
00200                    , G4RandGauss::shoot() * std::sqrt( kT*mass ) 
00201                    , G4RandGauss::shoot() * std::sqrt( kT*mass ) );
00202 
00203    theTarget.SetMomentum( v );
00204 
00205 
00206      G4ThreeVector the3Neutron = theNeutron.GetMomentum();
00207      G4double nEnergy = theNeutron.GetTotalEnergy();
00208      G4ThreeVector the3Target = theTarget.GetMomentum();
00209      G4double tEnergy = theTarget.GetTotalEnergy();
00210      G4ReactionProduct theCMS;
00211      G4double totE = nEnergy+tEnergy;
00212      G4ThreeVector the3CMS = the3Target+the3Neutron;
00213      theCMS.SetMomentum(the3CMS);
00214      G4double cmsMom = std::sqrt(the3CMS*the3CMS);
00215      G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
00216      theCMS.SetMass(sqrts);
00217      theCMS.SetTotalEnergy(totE);
00218 
00219        theNeutron.Lorentz(theNeutron, theCMS);
00220        theTarget.Lorentz(theTarget, theCMS);
00221        G4double en = theNeutron.GetTotalMomentum(); // already in CMS.
00222        G4ThreeVector cms3Mom=theNeutron.GetMomentum(); // for neutron direction in CMS
00223        G4double cms_theta=cms3Mom.theta();
00224        G4double cms_phi=cms3Mom.phi();
00225        G4ThreeVector tempVector;
00226        tempVector.setX(std::cos(theta)*std::sin(cms_theta)*std::cos(cms_phi)
00227                        +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::cos(cms_phi)
00228                        -std::sin(theta)*std::sin(phi)*std::sin(cms_phi)  );
00229        tempVector.setY(std::cos(theta)*std::sin(cms_theta)*std::sin(cms_phi)
00230                        +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::sin(cms_phi)
00231                        +std::sin(theta)*std::sin(phi)*std::cos(cms_phi)  );
00232        tempVector.setZ(std::cos(theta)*std::cos(cms_theta)
00233                        -std::sin(theta)*std::cos(phi)*std::sin(cms_theta)  );
00234        tempVector *= en;
00235        theNeutron.SetMomentum(tempVector);
00236        theTarget.SetMomentum(-tempVector);
00237        G4double tP = theTarget.GetTotalMomentum();
00238        G4double tM = theTarget.GetMass();
00239        theTarget.SetTotalEnergy(std::sqrt((tP+tM)*(tP+tM)-2.*tP*tM));
00240        theNeutron.Lorentz(theNeutron, -1.*theCMS);
00241        theTarget.Lorentz(theTarget, -1.*theCMS);
00242 
00243      theResult->SetEnergyChange(theNeutron.GetKineticEnergy());
00244      theResult->SetMomentumChange(theNeutron.GetMomentum().unit());
00245      G4DynamicParticle* theRecoil = new G4DynamicParticle;
00246 
00247      theRecoil->SetDefinition( G4ParticleTable::GetParticleTable()->FindIon( iZ, iA , 0, iZ ) );
00248      theRecoil->SetMomentum( theTarget.GetMomentum() );
00249 
00250      theResult->AddSecondary( theRecoil );
00251 
00252    return theResult; 
00253 
00254 }

Generated on Mon May 27 17:48:45 2013 for Geant4 by  doxygen 1.4.7