source: trunk/source/processes/hadronic/stopping/src/G4MuonMinusCaptureAtRest.cc@ 819

Last change on this file since 819 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 12.5 KB
Line 
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
58G4MuonMinusCaptureAtRest::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
72G4MuonMinusCaptureAtRest::~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
83G4bool G4MuonMinusCaptureAtRest::IsApplicable(const G4ParticleDefinition& p)
84{
85 return ( &p == G4MuonMinus::MuonMinus() );
86}
87
88//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89
90G4VParticleChange* 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
204G4ReactionProductVector* 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
Note: See TracBrowser for help on using the repository browser.