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
00038
00039
00040
00041
00042
00043
00044
00045 #include "G4MuMinusCapturePrecompound.hh"
00046 #include "Randomize.hh"
00047 #include "G4RandomDirection.hh"
00048 #include "G4PhysicalConstants.hh"
00049 #include "G4SystemOfUnits.hh"
00050 #include "G4MuonMinus.hh"
00051 #include "G4NeutrinoMu.hh"
00052 #include "G4Neutron.hh"
00053 #include "G4Proton.hh"
00054 #include "G4Triton.hh"
00055 #include "G4LorentzVector.hh"
00056 #include "G4ParticleDefinition.hh"
00057 #include "G4NucleiProperties.hh"
00058 #include "G4VPreCompoundModel.hh"
00059 #include "G4PreCompoundModel.hh"
00060 #include "G4HadronicInteractionRegistry.hh"
00061
00062
00063
00064 G4MuMinusCapturePrecompound::G4MuMinusCapturePrecompound(
00065 G4VPreCompoundModel* ptr)
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 }
00081
00082
00083
00084 G4MuMinusCapturePrecompound::~G4MuMinusCapturePrecompound()
00085 {
00086 result.Clear();
00087 }
00088
00089
00090
00091 G4HadFinalState*
00092 G4MuMinusCapturePrecompound::ApplyYourself(const G4HadProjectile& projectile,
00093 G4Nucleus& targetNucleus)
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
00108
00109
00110
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
00120
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
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
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154 G4LorentzVector momInitial(0.0,0.0,0.0,availableEnergy);
00155 G4LorentzVector momResidual, momNu;
00156
00157
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
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
00185 momNu.set(Enu*G4RandomDirection(), Enu);
00186
00187
00188 momNu.boost(bst);
00189 momResidual = momInitial - momNu;
00190 eEx = momResidual.mag() - residualMass;
00191
00192
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
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
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 }
00251
00252
00253
00254 void G4MuMinusCapturePrecompound::ModelDescription(std::ostream& outFile) const
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 }
00262
00263