#include <G4MuMinusCapturePrecompound.hh>
Inheritance diagram for G4MuMinusCapturePrecompound:
Public Member Functions | |
G4MuMinusCapturePrecompound (G4VPreCompoundModel *ptr=0) | |
~G4MuMinusCapturePrecompound () | |
G4HadFinalState * | ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus) |
void | ModelDescription (std::ostream &outFile) const |
Definition at line 65 of file G4MuMinusCapturePrecompound.hh.
G4MuMinusCapturePrecompound::G4MuMinusCapturePrecompound | ( | G4VPreCompoundModel * | ptr = 0 |
) |
Definition at line 64 of file G4MuMinusCapturePrecompound.cc.
References G4HadronicInteractionRegistry::FindModel(), G4ParticleDefinition::GetPDGMass(), G4HadronicInteractionRegistry::Instance(), G4MuonMinus::MuonMinus(), G4Neutron::Neutron(), and G4Proton::Proton().
00066 : G4HadronicInteraction("muMinusNuclearCapture") 00067 { 00068 fMuMass = G4MuonMinus::MuonMinus()->GetPDGMass(); 00069 fProton = G4Proton::Proton(); 00070 fNeutron = G4Neutron::Neutron(); 00071 fThreshold = 10*MeV; 00072 fPreCompound = ptr; 00073 if(!ptr) { 00074 G4HadronicInteraction* p = 00075 G4HadronicInteractionRegistry::Instance()->FindModel("PRECO"); 00076 ptr = static_cast<G4VPreCompoundModel*>(p); 00077 fPreCompound = ptr; 00078 if(!ptr) { fPreCompound = new G4PreCompoundModel(); } 00079 } 00080 }
G4MuMinusCapturePrecompound::~G4MuMinusCapturePrecompound | ( | ) |
Definition at line 84 of file G4MuMinusCapturePrecompound.cc.
References G4HadFinalState::Clear().
00085 { 00086 result.Clear(); 00087 }
G4HadFinalState * G4MuMinusCapturePrecompound::ApplyYourself | ( | const G4HadProjectile & | aTrack, | |
G4Nucleus & | targetNucleus | |||
) | [virtual] |
Implements G4HadronicInteraction.
Definition at line 92 of file G4MuMinusCapturePrecompound.cc.
References G4HadFinalState::Clear(), G4VPreCompoundModel::DeExcite(), FatalException, G4cout, G4endl, G4Exception(), G4RandomDirection(), G4UniformRand, G4DynamicParticle::Get4Momentum(), G4Nucleus::GetA_asInt(), G4HadProjectile::GetBoundEnergy(), G4HadProjectile::GetGlobalTime(), G4HadronicInteraction::GetModelName(), G4NucleiProperties::GetNuclearMass(), G4Fancy3DNucleus::GetNucleons(), G4HadFinalState::GetNumberOfSecondaries(), G4HadSecondary::GetParticle(), G4HadFinalState::GetSecondary(), G4Nucleus::GetZ_asInt(), G4Fancy3DNucleus::Init(), CLHEP::detail::n, G4NeutrinoMu::NeutrinoMu(), G4HadFinalState::SetStatusChange(), stopAndKill, G4Triton::Triton(), and G4HadronicInteraction::verboseLevel.
00094 { 00095 result.Clear(); 00096 result.SetStatusChange(stopAndKill); 00097 fTime = projectile.GetGlobalTime(); 00098 G4double time0 = fTime; 00099 00100 G4double muBindingEnergy = projectile.GetBoundEnergy(); 00101 00102 G4int Z = targetNucleus.GetZ_asInt(); 00103 G4int A = targetNucleus.GetA_asInt(); 00104 G4double massA = G4NucleiProperties::GetNuclearMass(A, Z); 00105 00106 /* 00107 G4cout << "G4MuMinusCapturePrecompound::ApplyYourself: Emu= " 00108 << muBindingEnergy << G4endl; 00109 */ 00110 // Energy on K-shell 00111 G4double muEnergy = fMuMass + muBindingEnergy; 00112 G4double muMom = std::sqrt(muBindingEnergy*(muBindingEnergy + 2.0*fMuMass)); 00113 G4double availableEnergy = massA + fMuMass - muBindingEnergy; 00114 G4double residualMass = G4NucleiProperties::GetNuclearMass(A, Z - 1); 00115 00116 G4ThreeVector vmu = muMom*G4RandomDirection(); 00117 G4LorentzVector aMuMom(vmu, muEnergy); 00118 00119 // p or 3He as a target 00120 // two body reaction mu- + A(Z,A) -> nuMu + A(Z-1,A) 00121 if((1 == Z && 1 == A) || (2 == Z && 3 == A)) { 00122 00123 G4ParticleDefinition* pd = 0; 00124 if(1 == Z) { pd = fNeutron; } 00125 else { pd = G4Triton::Triton(); } 00126 00127 // 00128 // Computation in assumption of CM reaction 00129 // 00130 G4double e = 0.5*(availableEnergy - 00131 residualMass*residualMass/availableEnergy); 00132 00133 G4ThreeVector nudir = G4RandomDirection(); 00134 AddNewParticle(G4NeutrinoMu::NeutrinoMu(), nudir, e); 00135 nudir *= -1.0; 00136 AddNewParticle(pd, nudir, availableEnergy - e - residualMass); 00137 00138 00139 } else { 00140 // sample mu- + p -> nuMu + n reaction in CM of muonic atom 00141 00142 // muon 00143 // 00144 // NOTE by K.Genser and J.Yarba: 00145 // The code below isn't working because emu always turns smaller than fMuMass 00146 // For this reason the sqrt is producing a NaN 00147 // 00148 // G4double emu = (availableEnergy*availableEnergy - massA*massA 00149 // + fMuMass*fMuMass)/(2*availableEnergy); 00150 // G4ThreeVector mudir = G4RandomDirection(); 00151 // G4LorentzVector momMuon(std::sqrt(emu*emu - fMuMass*fMuMass)*mudir, emu); 00152 00153 // nucleus 00154 G4LorentzVector momInitial(0.0,0.0,0.0,availableEnergy); 00155 G4LorentzVector momResidual, momNu; 00156 00157 // pick random proton inside nucleus 00158 G4double eEx; 00159 fNucleus.Init(A, Z); 00160 const std::vector<G4Nucleon>& nucleons= fNucleus.GetNucleons(); 00161 G4ParticleDefinition* pDef; 00162 00163 G4int nneutrons = 1; 00164 G4int reentryCount = 0; 00165 00166 do { 00167 ++reentryCount; 00168 G4int index = 0; 00169 do { 00170 index=G4int(A*G4UniformRand()); 00171 pDef = nucleons[index].GetDefinition(); 00172 } while(pDef != fProton); 00173 G4LorentzVector momP = nucleons[index].Get4Momentum(); 00174 00175 // Get CMS kinematics 00176 G4LorentzVector theCMS = momP + aMuMom; 00177 G4ThreeVector bst = theCMS.boostVector(); 00178 00179 G4double Ecms = theCMS.mag(); 00180 G4double Enu = 0.5*(Ecms - neutron_mass_c2*neutron_mass_c2/Ecms); 00181 eEx = 0.0; 00182 00183 if(Enu > 0.0) { 00184 // make the nu, and transform to lab; 00185 momNu.set(Enu*G4RandomDirection(), Enu); 00186 00187 // nu in lab. 00188 momNu.boost(bst); 00189 momResidual = momInitial - momNu; 00190 eEx = momResidual.mag() - residualMass; 00191 00192 // release neutron 00193 00194 if(eEx > 0.0) { 00195 G4double eth = residualMass - massA + fThreshold + 2*neutron_mass_c2; 00196 if(Ecms - Enu > eth) { 00197 theCMS -= momNu; 00198 G4double ekin = theCMS.e() - eth; 00199 G4ThreeVector dir = theCMS.vect().unit(); 00200 AddNewParticle(fNeutron, dir, ekin); 00201 momResidual -= 00202 result.GetSecondary(0)->GetParticle()->Get4Momentum(); 00203 --Z; 00204 --A; 00205 residualMass = G4NucleiProperties::GetNuclearMass(A, Z); 00206 nneutrons = 0; 00207 } 00208 } 00209 } 00210 if(Enu <= 0.0 && eEx <= 0.0 && reentryCount > 100) { 00211 G4ExceptionDescription ed; 00212 ed << "Call for " << GetModelName() << G4endl; 00213 ed << "Target Z= " << Z 00214 << " A= " << A << G4endl; 00215 ed << " ApplyYourself does not completed after 100 attempts" << G4endl; 00216 G4Exception("G4MuMinusCapturePrecompound::AtRestDoIt", "had006", 00217 FatalException, ed); 00218 } 00219 } while(eEx <= 0.0); 00220 00221 G4ThreeVector dir = momNu.vect().unit(); 00222 AddNewParticle(G4NeutrinoMu::NeutrinoMu(), dir, momNu.e()); 00223 00224 G4Fragment initialState(A, Z, momResidual); 00225 initialState.SetNumberOfExcitedParticle(nneutrons,0); 00226 initialState.SetNumberOfHoles(1,1); 00227 00228 // decay time for pre-compound/de-excitation starts from zero 00229 G4ReactionProductVector* rpv = fPreCompound->DeExcite(initialState); 00230 size_t n = rpv->size(); 00231 for(size_t i=0; i<n; ++i) { 00232 G4ReactionProduct* rp = (*rpv)[i]; 00233 00234 // reaction time 00235 fTime = time0 + rp->GetTOF(); 00236 G4ThreeVector direction = rp->GetMomentum().unit(); 00237 AddNewParticle(rp->GetDefinition(), direction, rp->GetKineticEnergy()); 00238 delete rp; 00239 } 00240 delete rpv; 00241 } 00242 if(verboseLevel > 1) 00243 G4cout << "G4MuMinusCapturePrecompound::ApplyYourself: Nsec= " 00244 << result.GetNumberOfSecondaries() 00245 <<" E0(MeV)= " <<availableEnergy/MeV 00246 <<" Mres(GeV)= " <<residualMass/GeV 00247 <<G4endl; 00248 00249 return &result; 00250 }
void G4MuMinusCapturePrecompound::ModelDescription | ( | std::ostream & | outFile | ) | const [virtual] |
Reimplemented from G4HadronicInteraction.
Definition at line 254 of file G4MuMinusCapturePrecompound.cc.
00255 { 00256 outFile << "Sampling of mu- capture by atomic nucleus from K-shell" 00257 << " mesoatom orbit.\n" 00258 << "Primary reaction mu- + p -> n + neutrino, neutron providing\n" 00259 << " initial excitation of the target nucleus and PreCompound" 00260 << " model samples final state\n"; 00261 }