Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 //-------------------------------------------- 28 // 29 // GEANT4 Class file 30 // 31 // File name: G4MuMinusCapturePrecompound 32 // 33 // Author: V.Ivanchenko (Vladimir.Ivant 34 // 35 // Creation date: 22 April 2012 on base of G4M 36 // 37 // 38 //-------------------------------------------- 39 // 40 // Modifications: 41 // 42 //-------------------------------------------- 43 44 #include "G4MuMinusCapturePrecompound.hh" 45 #include "Randomize.hh" 46 #include "G4RandomDirection.hh" 47 #include "G4PhysicalConstants.hh" 48 #include "G4SystemOfUnits.hh" 49 #include "G4MuonMinus.hh" 50 #include "G4NeutrinoMu.hh" 51 #include "G4Neutron.hh" 52 #include "G4Proton.hh" 53 #include "G4Triton.hh" 54 #include "G4LorentzVector.hh" 55 #include "G4ParticleDefinition.hh" 56 #include "G4NucleiProperties.hh" 57 #include "G4VPreCompoundModel.hh" 58 #include "G4PreCompoundModel.hh" 59 #include "G4HadronicInteractionRegistry.hh" 60 61 //....oooOO0OOooo........oooOO0OOooo........oo 62 63 G4MuMinusCapturePrecompound::G4MuMinusCaptureP 64 G4VPreCompoundModel* ptr) 65 : G4HadronicInteraction("muMinusNuclearCaptu 66 { 67 fMuMass = G4MuonMinus::MuonMinus()->GetPDGMa 68 fProton = G4Proton::Proton(); 69 fNeutron = G4Neutron::Neutron(); 70 fThreshold = 10*MeV; 71 fTime = 0.0; 72 fPreCompound = ptr; 73 if(!ptr) { 74 G4HadronicInteraction* p = 75 G4HadronicInteractionRegistry::Instance( 76 ptr = static_cast<G4VPreCompoundModel*>(p) 77 fPreCompound = ptr; 78 if(!ptr) { fPreCompound = new G4PreCompoun 79 } 80 } 81 82 //....oooOO0OOooo........oooOO0OOooo........oo 83 84 G4MuMinusCapturePrecompound::~G4MuMinusCapture 85 { 86 result.Clear(); 87 } 88 89 //....oooOO0OOooo........oooOO0OOooo........oo 90 91 G4HadFinalState* 92 G4MuMinusCapturePrecompound::ApplyYourself(con 93 G4Nucleus& targetNucleus) 94 { 95 result.Clear(); 96 result.SetStatusChange(stopAndKill); 97 fTime = projectile.GetGlobalTime(); 98 G4double time0 = fTime; 99 100 G4double muBindingEnergy = projectile.GetBou 101 102 G4int Z = targetNucleus.GetZ_asInt(); 103 G4int A = targetNucleus.GetA_asInt(); 104 G4double massA = G4NucleiProperties::GetNucl 105 106 /* 107 G4cout << "G4MuMinusCapturePrecompound::Appl 108 << muBindingEnergy << G4endl; 109 */ 110 // Energy on K-shell 111 G4double muEnergy = fMuMass + muBindingEnerg 112 G4double muMom =std::sqrt(muBindingEnergy*(m 113 G4double availableEnergy = massA + fMuMass - 114 G4double residualMass = G4NucleiProperties:: 115 116 G4ThreeVector vmu = muMom*G4RandomDirection( 117 G4LorentzVector aMuMom(vmu, muEnergy); 118 119 const G4double nenergy = keV; 120 121 // p or 3He as a target 122 // two body reaction mu- + A(Z,A) -> nuMu + 123 if((1 == Z && 1 == A) || (2 == Z && 3 == A)) 124 125 const G4ParticleDefinition* pd = 0; 126 if(1 == Z) { pd = fNeutron; } 127 else { pd = G4Triton::Triton(); } 128 129 // 130 // Computation in assumption of CM reacti 131 // 132 G4double e = 0.5*(availableEnergy - 133 residualMass*residualMass/availableE 134 135 G4ThreeVector nudir = G4RandomDirection(); 136 AddNewParticle(G4NeutrinoMu::NeutrinoMu(), 137 nudir *= -1.0; 138 AddNewParticle(pd, nudir, availableEnergy 139 140 // d or 4He as a target 141 // three body reaction mu- + A(Z,A) -> nuMu 142 // extra neutron produced at rest 143 } else if((1 == Z && 2 == A) || (2 == Z && 4 144 145 const G4ParticleDefinition* pd = 0; 146 if(1 == Z) { pd = fNeutron; } 147 else { pd = G4Triton::Triton(); } 148 149 availableEnergy -= neutron_mass_c2 - nener 150 residualMass = pd->GetPDGMass(); 151 152 // 153 // Computation in assumption of CM reacti 154 // 155 G4double e = 0.5*(availableEnergy - 156 residualMass*residualMass/availableE 157 158 G4ThreeVector nudir = G4RandomDirection(); 159 AddNewParticle(G4NeutrinoMu::NeutrinoMu(), 160 nudir *= -1.0; 161 AddNewParticle(pd, nudir, availableEnergy 162 163 // extra low-energy neutron 164 nudir = G4RandomDirection(); 165 AddNewParticle(fNeutron, nudir, nenergy); 166 167 } else { 168 // sample mu- + p -> nuMu + n reaction in 169 170 // nucleus 171 G4LorentzVector momInitial(0.0,0.0,0.0,ava 172 G4LorentzVector momResidual, momNu; 173 174 // pick random proton inside nucleus 175 G4double eEx; 176 fNucleus.Init(A, Z); 177 const std::vector<G4Nucleon>& nucleons= fN 178 const G4ParticleDefinition* pDef; 179 180 G4int reentryCount = 0; 181 182 do { 183 ++reentryCount; 184 G4int index = 0; 185 do { 186 index=G4int(A*G4UniformRand()); 187 pDef = nucleons[index].GetDefinition(); 188 } while(pDef != fProton); 189 G4LorentzVector momP = nucleons[index].G 190 191 // Get CMS kinematics 192 G4LorentzVector theCMS = momP + aMuMom; 193 G4ThreeVector bst = theCMS.boostVector() 194 195 G4double Ecms = theCMS.mag(); 196 G4double Enu = 0.5*(Ecms - neutron_mass 197 eEx = 0.0; 198 199 if(Enu > 0.0) { 200 // make the nu, and transform to lab; 201 momNu.set(Enu*G4RandomDirection(), Enu); 202 203 // nu in lab. 204 momNu.boost(bst); 205 momResidual = momInitial - momNu; 206 eEx = momResidual.mag() - residualMass; 207 if(eEx < 0.0 && eEx + nenergy >= 0.0) 208 momResidual.set(0.0, 0.0, 0.0, resid 209 eEx = 0.0; 210 } 211 } 212 // in the case of many iterations stop t 213 // with zero excitation energy 214 if(reentryCount > 100 && eEx < 0.0) { 215 G4ExceptionDescription ed; 216 ed << "Call for " << GetModelName() << G4end 217 ed << "Target Z= " << Z 218 << " A= " << A << " Eex(MeV)= " << eEx/ 219 ed << " ApplyYourself does not completed aft 220 << " excitation energy is set to zero"; 221 G4Exception("G4MuMinusCapturePrecompound::Ap 222 JustWarning, ed); 223 momResidual.set(0.0, 0.0, 0.0, residualMass) 224 eEx = 0.0; 225 } 226 // Loop checking, 06-Aug-2015, Vladimir 227 } while(eEx <= 0.0); 228 229 G4ThreeVector dir = momNu.vect().unit(); 230 AddNewParticle(G4NeutrinoMu::NeutrinoMu(), 231 232 G4Fragment initialState(A, Z-1, momResidua 233 initialState.SetNumberOfExcitedParticle(2, 234 initialState.SetNumberOfHoles(1,1); 235 236 // decay time for pre-compound/de-excitati 237 G4ReactionProductVector* rpv = fPreCompoun 238 size_t n = rpv->size(); 239 for(size_t i=0; i<n; ++i) { 240 G4ReactionProduct* rp = (*rpv)[i]; 241 242 // reaction time 243 fTime = time0 + rp->GetTOF(); 244 G4ThreeVector direction = rp->GetMomentu 245 AddNewParticle(rp->GetDefinition(), dire 246 delete rp; 247 } 248 delete rpv; 249 } 250 if(verboseLevel > 1) 251 G4cout << "G4MuMinusCapturePrecompound::Ap 252 << result.GetNumberOfSecondaries() 253 <<" E0(MeV)= " <<availableEnergy/MeV 254 <<" Mres(GeV)= " <<residualMass/GeV 255 <<G4endl; 256 257 return &result; 258 } 259 260 //....oooOO0OOooo........oooOO0OOooo........oo 261 262 void G4MuMinusCapturePrecompound::ModelDescrip 263 { 264 outFile << "Sampling of mu- capture by atomi 265 << " mesoatom orbit.\n" 266 << "Primary reaction mu- + p -> n + neutri 267 << " initial excitation of the target nuc 268 << " model samples final state\n"; 269 } 270 271 //....oooOO0OOooo........oooOO0OOooo........oo 272