[819] | 1 | // |
---|
| 2 | // ******************************************************************** |
---|
| 3 | // * License and Disclaimer * |
---|
| 4 | // * * |
---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
| 7 | // * conditions of the Geant4 Software License, included in the file * |
---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
| 9 | // * include a list of copyright holders. * |
---|
| 10 | // * * |
---|
| 11 | // * Neither the authors of this software system, nor their employing * |
---|
| 12 | // * institutes,nor the agencies providing financial support for this * |
---|
| 13 | // * work make any representation or warranty, express or implied, * |
---|
| 14 | // * regarding this software system or assume any liability for its * |
---|
| 15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
| 16 | // * for the full disclaimer and the limitation of liability. * |
---|
| 17 | // * * |
---|
| 18 | // * This code implementation is the result of the scientific and * |
---|
| 19 | // * technical work of the GEANT4 collaboration. * |
---|
| 20 | // * By using, copying, modifying or distributing the software (or * |
---|
| 21 | // * any work based on the software) you agree to acknowledge its * |
---|
| 22 | // * use in resulting scientific publications, and indicate your * |
---|
| 23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
| 24 | // ******************************************************************** |
---|
| 25 | // |
---|
[1055] | 26 | // $Id: G4MuonMinusCaptureAtRest.cc,v 1.54 2009/01/24 11:55:38 vnivanch Exp $ |
---|
| 27 | // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ |
---|
[819] | 28 | // |
---|
| 29 | // G4MuonMinusCaptureAtRest physics process |
---|
| 30 | // Larry Felawka (TRIUMF) and Art Olin (TRIUMF) |
---|
| 31 | // April 1998 |
---|
| 32 | //--------------------------------------------------------------------- |
---|
| 33 | // |
---|
| 34 | // Modifications: |
---|
[1055] | 35 | // 18/08/2000 V.Ivanchenko Update description |
---|
[819] | 36 | // 12/12/2003 H.P.Wellisch Completly rewrite mu-nuclear part |
---|
| 37 | // 17/05/2006 V.Ivanchenko Cleanup |
---|
| 38 | // 15/11/2006 V.Ivanchenko Review and rewrite all kinematics |
---|
| 39 | // 24/01/2007 V.Ivanchenko Force to work with integer Z and A |
---|
[1055] | 40 | // 23/01/2009 V.Ivanchenko Add deregistration |
---|
[819] | 41 | // |
---|
| 42 | //----------------------------------------------------------------------------- |
---|
| 43 | |
---|
| 44 | #include "G4MuonMinusCaptureAtRest.hh" |
---|
| 45 | #include "G4DynamicParticle.hh" |
---|
| 46 | #include "Randomize.hh" |
---|
| 47 | #include "G4He3.hh" |
---|
| 48 | #include "G4NeutrinoMu.hh" |
---|
| 49 | #include "G4Fragment.hh" |
---|
| 50 | #include "G4ReactionProductVector.hh" |
---|
| 51 | #include "G4Proton.hh" |
---|
| 52 | #include "G4PionPlus.hh" |
---|
| 53 | #include "G4GHEKinematicsVector.hh" |
---|
| 54 | #include "G4Fancy3DNucleus.hh" |
---|
| 55 | #include "G4ExcitationHandler.hh" |
---|
[962] | 56 | #include "G4HadronicProcessStore.hh" |
---|
[819] | 57 | |
---|
| 58 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 59 | |
---|
| 60 | G4MuonMinusCaptureAtRest::G4MuonMinusCaptureAtRest(const G4String& processName, |
---|
| 61 | G4ProcessType aType ) : |
---|
| 62 | G4VRestProcess (processName, aType), nCascade(0), targetZ(0), targetA(0), |
---|
| 63 | isInitialised(false) |
---|
| 64 | { |
---|
[962] | 65 | SetProcessSubType(fHadronAtRest); |
---|
[819] | 66 | Cascade = new G4GHEKinematicsVector [17]; |
---|
| 67 | pSelector = new G4StopElementSelector(); |
---|
| 68 | pEMCascade = new G4MuMinusCaptureCascade(); |
---|
| 69 | theN = new G4Fancy3DNucleus(); |
---|
| 70 | theHandler = new G4ExcitationHandler(); |
---|
[962] | 71 | G4HadronicProcessStore::Instance()->RegisterExtraProcess(this); |
---|
[819] | 72 | } |
---|
| 73 | |
---|
| 74 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 75 | |
---|
| 76 | G4MuonMinusCaptureAtRest::~G4MuonMinusCaptureAtRest() |
---|
| 77 | { |
---|
[1055] | 78 | G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this); |
---|
[819] | 79 | delete [] Cascade; |
---|
| 80 | delete pSelector; |
---|
| 81 | delete pEMCascade; |
---|
| 82 | delete theN; |
---|
| 83 | delete theHandler; |
---|
| 84 | } |
---|
| 85 | |
---|
| 86 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 87 | |
---|
| 88 | G4bool G4MuonMinusCaptureAtRest::IsApplicable(const G4ParticleDefinition& p) |
---|
| 89 | { |
---|
| 90 | return ( &p == G4MuonMinus::MuonMinus() ); |
---|
| 91 | } |
---|
| 92 | |
---|
| 93 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 94 | |
---|
[962] | 95 | void G4MuonMinusCaptureAtRest::PreparePhysicsTable(const G4ParticleDefinition& p) |
---|
| 96 | { |
---|
| 97 | G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this, &p); |
---|
| 98 | } |
---|
| 99 | |
---|
| 100 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 101 | |
---|
| 102 | void G4MuonMinusCaptureAtRest::BuildPhysicsTable(const G4ParticleDefinition& p) |
---|
| 103 | { |
---|
| 104 | G4HadronicProcessStore::Instance()->PrintInfo(&p); |
---|
| 105 | } |
---|
| 106 | |
---|
| 107 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 108 | |
---|
[819] | 109 | G4VParticleChange* G4MuonMinusCaptureAtRest::AtRestDoIt(const G4Track& track, |
---|
| 110 | const G4Step&) |
---|
| 111 | { |
---|
| 112 | // |
---|
| 113 | // Handles MuonMinuss at rest; a MuonMinus can either create secondaries or |
---|
| 114 | // do nothing (in which case it should be sent back to decay-handling |
---|
| 115 | // section |
---|
| 116 | // |
---|
| 117 | aParticleChange.Initialize(track); |
---|
| 118 | |
---|
| 119 | // select element and get Z,A. |
---|
| 120 | G4Element* aEle = pSelector->GetElement(track.GetMaterial()); |
---|
| 121 | targetZ = aEle->GetZ(); |
---|
| 122 | targetA = G4double(G4int(aEle->GetN()+0.5)); |
---|
| 123 | G4int ni = 0; |
---|
| 124 | |
---|
| 125 | G4IsotopeVector* isv = aEle->GetIsotopeVector(); |
---|
| 126 | if(isv) ni = isv->size(); |
---|
| 127 | |
---|
| 128 | if(ni == 1) { |
---|
| 129 | targetA = G4double(aEle->GetIsotope(0)->GetN()); |
---|
| 130 | } else if(ni > 1) { |
---|
| 131 | G4double* ab = aEle->GetRelativeAbundanceVector(); |
---|
| 132 | G4double y = G4UniformRand(); |
---|
| 133 | G4int j = -1; |
---|
| 134 | ni--; |
---|
| 135 | do { |
---|
| 136 | j++; |
---|
| 137 | y -= ab[j]; |
---|
| 138 | } while (y > 0.0 && j < ni); |
---|
| 139 | targetA = G4double(aEle->GetIsotope(j)->GetN()); |
---|
| 140 | } |
---|
| 141 | |
---|
| 142 | // Do the electromagnetic cascade of the muon in the nuclear field. |
---|
| 143 | nCascade = 0; |
---|
| 144 | targetMass = G4NucleiProperties::GetNuclearMass(targetA, targetZ); |
---|
| 145 | nCascade = pEMCascade->DoCascade(targetZ, targetMass, Cascade); |
---|
| 146 | |
---|
| 147 | // Decide on Decay or Capture, and doit. |
---|
| 148 | G4double lambdac = pSelector->GetMuonCaptureRate(targetZ, targetA); |
---|
| 149 | G4double lambdad = pSelector->GetMuonDecayRate(targetZ, targetA); |
---|
| 150 | G4double lambda = lambdac + lambdad; |
---|
| 151 | |
---|
| 152 | // === Throw for capture time. |
---|
| 153 | |
---|
| 154 | G4double tDelay = -std::log(G4UniformRand()) / lambda; |
---|
| 155 | |
---|
| 156 | G4ReactionProductVector * captureResult = 0; |
---|
| 157 | G4int nEmSecondaries = nCascade; |
---|
| 158 | G4int nSecondaries = nCascade; |
---|
| 159 | /* |
---|
| 160 | G4cout << "lambda= " << lambda << " lambdac= " << lambdac |
---|
| 161 | << " nem= " << nEmSecondaries << G4endl; |
---|
| 162 | */ |
---|
| 163 | if( G4UniformRand()*lambda > lambdac) |
---|
| 164 | pEMCascade->DoBoundMuonMinusDecay(targetZ, &nEmSecondaries, Cascade); |
---|
| 165 | else |
---|
| 166 | captureResult = DoMuCapture(); |
---|
| 167 | |
---|
| 168 | // fill the final state |
---|
| 169 | if(captureResult) nSecondaries += captureResult->size(); |
---|
| 170 | else nSecondaries = nEmSecondaries; |
---|
| 171 | //G4cout << " nsec= " << nSecondaries << " nem= " << nEmSecondaries << G4endl; |
---|
| 172 | |
---|
| 173 | aParticleChange.SetNumberOfSecondaries( nSecondaries ); |
---|
| 174 | |
---|
| 175 | G4double globalTime = track.GetGlobalTime(); |
---|
| 176 | G4ThreeVector position = track.GetPosition(); |
---|
| 177 | // Store nuclear cascade |
---|
| 178 | if(captureResult) { |
---|
| 179 | G4int n = captureResult->size(); |
---|
| 180 | for ( G4int isec = 0; isec < n; isec++ ) { |
---|
| 181 | G4ReactionProduct* aParticle = captureResult->operator[](isec); |
---|
| 182 | G4DynamicParticle * aNewParticle = new G4DynamicParticle(); |
---|
| 183 | aNewParticle->SetDefinition( aParticle->GetDefinition() ); |
---|
| 184 | G4LorentzVector itV(aParticle->GetTotalEnergy(), aParticle->GetMomentum()); |
---|
| 185 | aNewParticle->SetMomentum(itV.vect()); |
---|
| 186 | G4double localtime = globalTime + tDelay + aParticle->GetTOF(); |
---|
| 187 | G4Track* aNewTrack = new G4Track( aNewParticle, localtime, position); |
---|
| 188 | aNewTrack->SetTouchableHandle(track.GetTouchableHandle()); |
---|
| 189 | aParticleChange.AddSecondary( aNewTrack ); |
---|
| 190 | delete aParticle; |
---|
| 191 | } |
---|
| 192 | delete captureResult; |
---|
| 193 | } |
---|
| 194 | |
---|
| 195 | // Store electromagnetic cascade |
---|
| 196 | |
---|
| 197 | if(nEmSecondaries > 0) { |
---|
| 198 | |
---|
| 199 | for ( G4int isec = 0; isec < nEmSecondaries; isec++ ) { |
---|
| 200 | G4ParticleDefinition* pd = Cascade[isec].GetParticleDef(); |
---|
| 201 | G4double localtime = globalTime; |
---|
| 202 | if(isec >= nCascade) localtime += tDelay; |
---|
| 203 | if(pd) { |
---|
| 204 | G4DynamicParticle* aNewParticle = new G4DynamicParticle; |
---|
| 205 | aNewParticle->SetDefinition( pd ); |
---|
| 206 | aNewParticle->SetMomentum( Cascade[isec].GetMomentum() ); |
---|
| 207 | |
---|
| 208 | G4Track* aNewTrack = new G4Track( aNewParticle, localtime, position ); |
---|
| 209 | aNewTrack->SetTouchableHandle(track.GetTouchableHandle()); |
---|
| 210 | aParticleChange.AddSecondary( aNewTrack ); |
---|
| 211 | } |
---|
| 212 | } |
---|
| 213 | } |
---|
| 214 | |
---|
| 215 | aParticleChange.ProposeLocalEnergyDeposit(0.0); |
---|
| 216 | aParticleChange.ProposeTrackStatus(fStopAndKill); |
---|
| 217 | |
---|
| 218 | return &aParticleChange; |
---|
| 219 | } |
---|
| 220 | |
---|
| 221 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 222 | |
---|
| 223 | G4ReactionProductVector* G4MuonMinusCaptureAtRest::DoMuCapture() |
---|
| 224 | { |
---|
| 225 | G4double mumass = G4MuonMinus::MuonMinus()->GetPDGMass(); |
---|
| 226 | G4double muBindingEnergy = pEMCascade->GetKShellEnergy(targetZ); |
---|
| 227 | /* |
---|
| 228 | G4cout << "G4MuonMinusCaptureAtRest::DoMuCapture called Emu= " |
---|
| 229 | << muBindingEnergy << G4endl; |
---|
| 230 | */ |
---|
| 231 | // Energy on K-shell |
---|
| 232 | G4double muEnergy = mumass + muBindingEnergy; |
---|
| 233 | G4double muMom = std::sqrt(muBindingEnergy*(muBindingEnergy + 2.0*mumass)); |
---|
| 234 | G4double availableEnergy = targetMass + mumass - muBindingEnergy; |
---|
| 235 | G4LorentzVector momInitial(0.0,0.0,0.0,availableEnergy); |
---|
| 236 | G4LorentzVector momResidual; |
---|
| 237 | |
---|
| 238 | G4ThreeVector vmu = muMom*pEMCascade->GetRandomVec(); |
---|
| 239 | G4LorentzVector aMuMom(vmu, muEnergy); |
---|
| 240 | |
---|
| 241 | G4double residualMass = |
---|
| 242 | G4NucleiProperties::GetNuclearMass(targetA, targetZ - 1.0); |
---|
| 243 | |
---|
| 244 | G4ReactionProductVector* aPreResult = 0; |
---|
| 245 | G4ReactionProduct* aNu = new G4ReactionProduct(); |
---|
| 246 | aNu->SetDefinition( G4NeutrinoMu::NeutrinoMu() ); |
---|
| 247 | |
---|
| 248 | G4int iz = G4int(targetZ); |
---|
| 249 | G4int ia = G4int(targetA); |
---|
| 250 | |
---|
| 251 | // proton as a target |
---|
| 252 | if(iz <= 2) { |
---|
| 253 | |
---|
| 254 | if(ia > 1) { |
---|
| 255 | if(iz == 1 && ia == 2) { |
---|
| 256 | availableEnergy -= neutron_mass_c2; |
---|
| 257 | } else if(iz == 1 && ia == 3) { |
---|
| 258 | availableEnergy -= 2.0*neutron_mass_c2; |
---|
| 259 | } else if(iz == 2) { |
---|
| 260 | G4ParticleDefinition* pd = 0; |
---|
| 261 | if(ia == 3) pd = G4Deuteron::Deuteron(); |
---|
| 262 | if(ia == 4) pd = G4Triton::Triton(); |
---|
| 263 | else |
---|
| 264 | pd = G4ParticleTable::GetParticleTable()->FindIon(1,ia-1,0,1); |
---|
| 265 | |
---|
| 266 | // G4cout << "Extra " << pd->GetParticleName() << G4endl; |
---|
| 267 | availableEnergy -= pd->GetPDGMass(); |
---|
| 268 | } |
---|
| 269 | } |
---|
| 270 | // |
---|
| 271 | // Computation in assumption of CM collision of mu and nucleaon |
---|
| 272 | // |
---|
| 273 | G4double Enu = 0.5*(availableEnergy - |
---|
| 274 | neutron_mass_c2*neutron_mass_c2/availableEnergy); |
---|
| 275 | |
---|
| 276 | // make the nu, and transform to lab; |
---|
| 277 | G4ThreeVector nu3Mom = Enu*pEMCascade->GetRandomVec(); |
---|
| 278 | |
---|
| 279 | G4ReactionProduct* aN = new G4ReactionProduct(); |
---|
| 280 | aN->SetDefinition( G4Neutron::Neutron() ); |
---|
| 281 | aN->SetTotalEnergy( availableEnergy - Enu ); |
---|
| 282 | aN->SetMomentum( -nu3Mom ); |
---|
| 283 | |
---|
| 284 | aNu->SetTotalEnergy( Enu ); |
---|
| 285 | aNu->SetMomentum( nu3Mom ); |
---|
| 286 | aPreResult = new G4ReactionProductVector(); |
---|
| 287 | |
---|
| 288 | aPreResult->push_back(aN ); |
---|
| 289 | aPreResult->push_back(aNu); |
---|
| 290 | |
---|
| 291 | if(verboseLevel > 1) |
---|
| 292 | G4cout << "DoMuCapture on H or He" |
---|
| 293 | <<" EkinN(MeV)= " << (availableEnergy - Enu - neutron_mass_c2)/MeV |
---|
| 294 | <<" Enu(MeV)= "<<aNu->GetTotalEnergy()/MeV |
---|
| 295 | <<" n= " << aPreResult->size() |
---|
| 296 | <<G4endl; |
---|
| 297 | |
---|
| 298 | return aPreResult; |
---|
| 299 | } |
---|
| 300 | |
---|
| 301 | // pick random proton inside nucleus |
---|
| 302 | G4double eEx; |
---|
| 303 | do { |
---|
| 304 | theN->Init(targetA, targetZ); |
---|
| 305 | G4LorentzVector thePMom; |
---|
| 306 | G4Nucleon * aNucleon = 0; |
---|
| 307 | G4int theProtonCounter = G4int( targetZ * G4UniformRand() ); |
---|
| 308 | G4int counter = 0; |
---|
| 309 | theN->StartLoop(); |
---|
| 310 | |
---|
| 311 | while( (aNucleon=theN->GetNextNucleon()) ) { |
---|
| 312 | |
---|
| 313 | if( aNucleon->GetDefinition() == G4Proton::Proton() ) { |
---|
| 314 | counter++; |
---|
| 315 | if(counter == theProtonCounter) { |
---|
| 316 | thePMom = aNucleon->GetMomentum(); |
---|
| 317 | break; |
---|
| 318 | } |
---|
| 319 | } |
---|
| 320 | } |
---|
| 321 | |
---|
| 322 | // Get the nu momentum in the CMS |
---|
| 323 | G4LorentzVector theCMS = thePMom + aMuMom; |
---|
| 324 | G4ThreeVector bst = theCMS.boostVector(); |
---|
| 325 | |
---|
| 326 | G4double Ecms = theCMS.mag(); |
---|
| 327 | G4double Enu = 0.5*(Ecms - neutron_mass_c2*neutron_mass_c2/Ecms); |
---|
| 328 | eEx = 0.0; |
---|
| 329 | |
---|
| 330 | if(Enu > 0.0) { |
---|
| 331 | // make the nu, and transform to lab; |
---|
| 332 | G4ThreeVector nu3Mom = Enu*pEMCascade->GetRandomVec(); |
---|
| 333 | G4LorentzVector nuMom(nu3Mom, Enu); |
---|
| 334 | |
---|
| 335 | // nu in lab. |
---|
| 336 | nuMom.boost(bst); |
---|
| 337 | aNu->SetTotalEnergy( nuMom.e() ); |
---|
| 338 | aNu->SetMomentum( nuMom.vect() ); |
---|
| 339 | |
---|
| 340 | // make residual |
---|
| 341 | momResidual = momInitial - nuMom; |
---|
| 342 | |
---|
| 343 | // Call pre-compound on the rest. |
---|
| 344 | eEx = momResidual.mag(); |
---|
| 345 | if(verboseLevel > 1) |
---|
| 346 | G4cout << "G4MuonMinusCaptureAtRest::DoMuCapture: " |
---|
| 347 | << " Eex(MeV)= " << (eEx-residualMass)/MeV |
---|
| 348 | << " Enu(MeV)= "<<aNu->GetTotalEnergy()/MeV |
---|
| 349 | <<G4endl; |
---|
| 350 | } |
---|
| 351 | } while(eEx <= residualMass); |
---|
| 352 | |
---|
| 353 | // G4cout << "muonCapture : " << eEx << " " << residualMass |
---|
| 354 | // << " A,Z= " << targetA << ", "<< targetZ |
---|
| 355 | // << " " << G4int(targetA) << ", " << G4int(targetZ) << G4endl; |
---|
| 356 | |
---|
| 357 | // |
---|
| 358 | // Start Deexcitation |
---|
| 359 | // |
---|
| 360 | G4ThreeVector fromBreit = momResidual.boostVector(); |
---|
| 361 | G4LorentzVector fscm(0.0,0.0,0.0, eEx); |
---|
| 362 | G4Fragment anInitialState; |
---|
| 363 | anInitialState.SetA(targetA); |
---|
| 364 | anInitialState.SetZ(G4double(iz - 1)); |
---|
| 365 | anInitialState.SetNumberOfParticles(2); |
---|
| 366 | anInitialState.SetNumberOfCharged(0); |
---|
| 367 | anInitialState.SetNumberOfHoles(1); |
---|
| 368 | anInitialState.SetMomentum(fscm); |
---|
| 369 | aPreResult = theHandler->BreakItUp(anInitialState); |
---|
| 370 | |
---|
| 371 | G4ReactionProductVector::iterator ires; |
---|
| 372 | G4double eBal = availableEnergy - aNu->GetTotalEnergy(); |
---|
| 373 | for(ires=aPreResult->begin(); ires!=aPreResult->end(); ires++) { |
---|
| 374 | G4LorentzVector itV((*ires)->GetTotalEnergy(), (*ires)->GetMomentum()); |
---|
| 375 | itV.boost(fromBreit); |
---|
| 376 | (*ires)->SetTotalEnergy(itV.t()); |
---|
| 377 | (*ires)->SetMomentum(itV.vect()); |
---|
| 378 | eBal -= itV.t(); |
---|
| 379 | } |
---|
| 380 | // |
---|
| 381 | // fill neutrino into result |
---|
| 382 | // |
---|
| 383 | aPreResult->push_back(aNu); |
---|
| 384 | |
---|
| 385 | if(verboseLevel > 1) |
---|
| 386 | G4cout << "DoMuCapture: Nsec= " |
---|
| 387 | << aPreResult->size() << " Ebalance(MeV)= " << eBal/MeV |
---|
| 388 | <<" E0(MeV)= " <<availableEnergy/MeV |
---|
| 389 | <<" Mres(GeV)= " <<residualMass/GeV |
---|
| 390 | <<G4endl; |
---|
| 391 | |
---|
| 392 | return aPreResult; |
---|
| 393 | } |
---|
| 394 | |
---|