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
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #include "G4NeutronRadCapture.hh"
00038 #include "G4SystemOfUnits.hh"
00039 #include "G4ParticleDefinition.hh"
00040 #include "G4Fragment.hh"
00041 #include "G4FragmentVector.hh"
00042 #include "G4NucleiProperties.hh"
00043 #include "G4PhotonEvaporation.hh"
00044 #include "G4DynamicParticle.hh"
00045 #include "G4ParticleDefinition.hh"
00046 #include "G4ParticleTable.hh"
00047 #include "G4IonTable.hh"
00048 #include "G4Deuteron.hh"
00049 #include "G4Triton.hh"
00050 #include "G4He3.hh"
00051 #include "G4Alpha.hh"
00052
00053 G4NeutronRadCapture::G4NeutronRadCapture()
00054 : G4HadronicInteraction("nRadCapture")
00055 {
00056 lowestEnergyLimit = 0.1*eV;
00057 SetMinEnergy( 0.0*GeV );
00058 SetMaxEnergy( 100.*TeV );
00059 photonEvaporation = new G4PhotonEvaporation();
00060
00061 }
00062
00063 G4NeutronRadCapture::~G4NeutronRadCapture()
00064 {
00065 delete photonEvaporation;
00066 }
00067
00068 G4HadFinalState* G4NeutronRadCapture::ApplyYourself(
00069 const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
00070 {
00071 theParticleChange.Clear();
00072 theParticleChange.SetStatusChange(stopAndKill);
00073
00074 G4int A = theNucleus.GetA_asInt();
00075 G4int Z = theNucleus.GetZ_asInt();
00076
00077
00078 G4double m1 = G4NucleiProperties::GetNuclearMass(A, Z);
00079 G4LorentzVector lv0(0.0,0.0,0.0,m1);
00080 G4LorentzVector lv1 = aTrack.Get4Momentum() + lv0;
00081
00082
00083 if(A <= 3) {
00084
00085 G4ThreeVector bst = lv1.boostVector();
00086 G4double M = lv1.mag();
00087
00088 ++A;
00089 G4double mass = G4NucleiProperties::GetNuclearMass(A, Z);
00090 if(M - mass <= lowestEnergyLimit) {
00091 return &theParticleChange;
00092 }
00093
00094 if (verboseLevel > 1) {
00095 G4cout << "G4NeutronRadCapture::DoIt: Eini(MeV)="
00096 << aTrack.GetKineticEnergy()/MeV << " Eexc(MeV)= "
00097 << (M - mass)/MeV
00098 << " Z= " << Z << " A= " << A << G4endl;
00099 }
00100 G4double e1 = (M - mass)*(M + mass)/(2*M);
00101 G4double cost = 2.0*G4UniformRand() - 1.0;
00102 if(cost > 1.0) {cost = 1.0;}
00103 else if(cost < -1.0) {cost = -1.0;}
00104 G4double sint = std::sqrt((1. - cost)*(1.0 + cost));
00105 G4double phi = G4UniformRand()*CLHEP::twopi;
00106 G4LorentzVector lv2(e1*sint*std::cos(phi),e1*sint*std::sin(phi),e1*cost,e1);
00107 lv2.boost(bst);
00108 theParticleChange.AddSecondary(new G4DynamicParticle(G4Gamma::Gamma(), lv2));
00109 G4ParticleDefinition* theDef = 0;
00110
00111 lv1 -= lv2;
00112 if (Z == 1 && A == 2) {theDef = G4Deuteron::Deuteron();}
00113 else if (Z == 1 && A == 3) {theDef = G4Triton::Triton();}
00114 else if (Z == 2 && A == 3) {theDef = G4He3::He3();}
00115 else if (Z == 2 && A == 4) {theDef = G4Alpha::Alpha();}
00116 else
00117 {
00118 theDef =
00119 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(Z,A,0.0);
00120 }
00121
00122 if (verboseLevel > 1) {
00123 G4cout << "Gamma 4-mom: " << lv2 << " "
00124 << theDef->GetParticleName() << " " << lv1 << G4endl;
00125 }
00126 if(theDef) {
00127 theParticleChange.AddSecondary(new G4DynamicParticle(theDef, lv1));
00128 }
00129
00130
00131 } else {
00132
00133 G4Fragment* aFragment = new G4Fragment(A+1, Z, lv1);
00134
00135 if (verboseLevel > 1) {
00136 G4cout << "G4NeutronRadCapture::ApplyYourself initial G4Fragmet:" << G4endl;
00137 G4cout << aFragment << G4endl;
00138 }
00139
00140
00141
00142
00143 G4FragmentVector* fv = photonEvaporation->BreakUpFragment(aFragment);
00144 if(!fv) { fv = new G4FragmentVector(); }
00145 fv->push_back(aFragment);
00146 size_t n = fv->size();
00147
00148 if (verboseLevel > 1) {
00149 G4cout << "G4NeutronRadCapture: " << n << " final particle" << G4endl;
00150 }
00151 for(size_t i=0; i<n; ++i) {
00152 G4Fragment* f = (*fv)[i];
00153 G4LorentzVector lvres = f->GetMomentum();
00154 Z = f->GetZ_asInt();
00155 A = f->GetA_asInt();
00156
00157 G4ParticleDefinition* theDef = 0;
00158 if(0 == Z && 0 == A) {theDef = f->GetParticleDefinition();}
00159 else
00160 {
00161 theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(Z,A,0.0);
00162 }
00163
00164 if (verboseLevel > 1) {
00165 G4cout << i << ". " << theDef->GetParticleName()
00166 << " " << lvres << G4endl;
00167 }
00168 if(theDef) {
00169 theParticleChange.AddSecondary(new G4DynamicParticle(theDef, lvres));
00170 }
00171 delete f;
00172 }
00173 delete fv;
00174 }
00175 return &theParticleChange;
00176 }
00177