G4LENDElastic.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 
00027 #include "G4LENDElastic.hh"
00028 #include "G4PhysicalConstants.hh"
00029 #include "G4SystemOfUnits.hh"
00030 #include "G4Nucleus.hh"
00031 #include "G4ParticleTable.hh"
00032 
00033 G4HadFinalState * G4LENDElastic::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& aTarg )
00034 {
00035 
00036    G4double temp = aTrack.GetMaterial()->GetTemperature();
00037 
00038    //G4int iZ = int ( aTarg.GetZ() );
00039    //G4int iA = int ( aTarg.GetN() );
00040    //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this) 
00041    G4int iZ = aTarg.GetZ_asInt();
00042    G4int iA = aTarg.GetA_asInt();
00043 
00044    G4double ke = aTrack.GetKineticEnergy();
00045 
00046    //G4HadFinalState* theResult = new G4HadFinalState();
00047    G4HadFinalState* theResult = &theParticleChange;
00048    theResult->Clear();
00049 
00050    G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA ) )->second->GetTarget();
00051    G4double aMu = aTarget->getElasticFinalState( ke*MeV, temp, NULL, NULL );
00052 
00053    G4double phi = twopi*G4UniformRand();
00054    G4double theta = std::acos( aMu );
00055    //G4double sinth = std::sin( theta );
00056 
00057    G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>( aTrack.GetDefinition() ) );
00058    theNeutron.SetMomentum( aTrack.Get4Momentum().vect() );
00059    theNeutron.SetKineticEnergy( ke );
00060 
00061 //G4cout << "iZ " << iZ << " iA " << iA  << G4endl;
00062 
00063    G4ReactionProduct theTarget( G4ParticleTable::GetParticleTable()->FindIon( iZ , iA , 0 , iZ ) );
00064 
00065    G4double mass = G4ParticleTable::GetParticleTable()->FindIon( iZ , iA , 0 , iZ )->GetPDGMass();
00066 
00067 // add Thermal motion 
00068    G4double kT = k_Boltzmann*temp;
00069    G4ThreeVector v ( G4RandGauss::shoot() * std::sqrt( kT*mass ) 
00070                    , G4RandGauss::shoot() * std::sqrt( kT*mass ) 
00071                    , G4RandGauss::shoot() * std::sqrt( kT*mass ) );
00072    theTarget.SetMomentum( v );
00073 
00074      G4ThreeVector the3Neutron = theNeutron.GetMomentum();
00075      G4double nEnergy = theNeutron.GetTotalEnergy();
00076      G4ThreeVector the3Target = theTarget.GetMomentum();
00077      G4double tEnergy = theTarget.GetTotalEnergy();
00078      G4ReactionProduct theCMS;
00079      G4double totE = nEnergy+tEnergy;
00080      G4ThreeVector the3CMS = the3Target+the3Neutron;
00081      theCMS.SetMomentum(the3CMS);
00082      G4double cmsMom = std::sqrt(the3CMS*the3CMS);
00083      G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
00084      theCMS.SetMass(sqrts);
00085      theCMS.SetTotalEnergy(totE);
00086 
00087        theNeutron.Lorentz(theNeutron, theCMS);
00088        theTarget.Lorentz(theTarget, theCMS);
00089        G4double en = theNeutron.GetTotalMomentum(); // already in CMS.
00090        G4ThreeVector cms3Mom=theNeutron.GetMomentum(); // for neutron direction in CMS
00091        G4double cms_theta=cms3Mom.theta();
00092        G4double cms_phi=cms3Mom.phi();
00093        G4ThreeVector tempVector;
00094        tempVector.setX( std::cos(theta)*std::sin(cms_theta)*std::cos(cms_phi)
00095                        +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::cos(cms_phi)
00096                        -std::sin(theta)*std::sin(phi)*std::sin(cms_phi) );
00097        tempVector.setY( std::cos(theta)*std::sin(cms_theta)*std::sin(cms_phi)
00098                        +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::sin(cms_phi)
00099                        +std::sin(theta)*std::sin(phi)*std::cos(cms_phi) );
00100        tempVector.setZ( std::cos(theta)*std::cos(cms_theta)
00101                        -std::sin(theta)*std::cos(phi)*std::sin(cms_theta) );
00102        tempVector *= en;
00103        theNeutron.SetMomentum(tempVector);
00104        theTarget.SetMomentum(-tempVector);
00105        G4double tP = theTarget.GetTotalMomentum();
00106        G4double tM = theTarget.GetMass();
00107        theTarget.SetTotalEnergy(std::sqrt((tP+tM)*(tP+tM)-2.*tP*tM));
00108 
00109 
00110        theNeutron.Lorentz(theNeutron, -1.*theCMS);
00111 
00112 //110913 Add Protection for very low energy (1e-6eV) scattering 
00113       if ( theNeutron.GetKineticEnergy() <= 0 )
00114       {
00115          theNeutron.SetTotalEnergy ( theNeutron.GetMass() * ( 1 + std::pow( 10 , -15.65 ) ) );
00116       }
00117 
00118       theTarget.Lorentz(theTarget, -1.*theCMS);
00119       if ( theTarget.GetKineticEnergy() < 0 )
00120       {
00121          theTarget.SetTotalEnergy ( theTarget.GetMass() * ( 1 + std::pow( 10 , -15.65 ) ) );
00122       }
00123 //110913 END
00124 
00125        theTarget.Lorentz(theTarget, -1.*theCMS);
00126 
00127      theResult->SetEnergyChange(theNeutron.GetKineticEnergy());
00128      theResult->SetMomentumChange(theNeutron.GetMomentum().unit());
00129      G4DynamicParticle* theRecoil = new G4DynamicParticle;
00130 
00131 //     theRecoil->SetDefinition( ionTable->GetIon( iZ , iA ) ); 
00132        theRecoil->SetDefinition( G4ParticleTable::GetParticleTable()->FindIon( iZ, iA , 0, iZ ));
00133      theRecoil->SetMomentum( theTarget.GetMomentum() );
00134 
00135      theResult->AddSecondary( theRecoil );
00136 
00137    return theResult; 
00138 
00139 }
00140 

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