00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "G4LENDInelastic.hh"
00027 #include "G4SystemOfUnits.hh"
00028 #include "G4Nucleus.hh"
00029 #include "G4ParticleTable.hh"
00030
00031 G4HadFinalState * G4LENDInelastic::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& aTarg )
00032 {
00033
00034 G4ThreeVector proj_p = aTrack.Get4Momentum().vect();
00035
00036 G4double temp = aTrack.GetMaterial()->GetTemperature();
00037
00038
00039
00040
00041 G4int iZ = aTarg.GetZ_asInt();
00042 G4int iA = aTarg.GetA_asInt();
00043
00044
00045 G4double ke = aTrack.GetKineticEnergy();
00046
00047
00048 G4HadFinalState* theResult = &theParticleChange;
00049 theResult->Clear();
00050
00051 G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA ) )->second->GetTarget();
00052 std::vector<G4GIDI_Product>* products = aTarget->getOthersFinalState( ke*MeV, temp, NULL, NULL );
00053 if ( products != NULL )
00054 {
00055
00056 G4ThreeVector psum(0);
00057
00058 int totN = 0;
00059 for ( G4int j = 0; j < int( products->size() ); j++ )
00060 {
00061
00062 G4int jZ = (*products)[j].Z;
00063 G4int jA = (*products)[j].A;
00064
00065
00066
00067
00068
00069
00070
00071
00072 G4DynamicParticle* theSec = new G4DynamicParticle;
00073
00074 if ( jA == 1 && jZ == 1 )
00075 {
00076 theSec->SetDefinition( G4Proton::Proton() );
00077 totN += 1;
00078 }
00079 else if ( jA == 1 && jZ == 0 )
00080 {
00081 theSec->SetDefinition( G4Neutron::Neutron() );
00082 totN += 1;
00083 }
00084 else if ( jZ > 0 )
00085 {
00086 if ( jA != 0 )
00087 {
00088 theSec->SetDefinition( G4ParticleTable::GetParticleTable()->FindIon( jZ , jA , 0 , 0 ) );
00089 totN += jA;
00090 }
00091 else
00092 {
00093 theSec->SetDefinition( G4ParticleTable::GetParticleTable()->FindIon( jZ , iA+1-totN , 0 , 0 ) );
00094 }
00095 }
00096 else
00097 {
00098 theSec->SetDefinition( G4Gamma::Gamma() );
00099 }
00100
00101 G4ThreeVector p( (*products)[j].px*MeV , (*products)[j].py*MeV , (*products)[j].pz*MeV );
00102 psum += p;
00103 if ( p.mag() == 0 ) p = proj_p - psum;
00104
00105 theSec->SetMomentum( p );
00106
00107 theResult->AddSecondary( theSec );
00108 }
00109 }
00110 delete products;
00111
00112 theResult->SetStatusChange( stopAndKill );
00113
00114 return theResult;
00115
00116 }