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