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

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

maj sur la beta de geant 4.9.3
File size: 13.2 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.54 2009/01/24 11:55:38 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
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// 23/01/2009 V.Ivanchenko Add deregistration
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"
56#include "G4HadronicProcessStore.hh"
57
58//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
59
60G4MuonMinusCaptureAtRest::G4MuonMinusCaptureAtRest(const G4String& processName,
61 G4ProcessType aType ) :
62 G4VRestProcess (processName, aType), nCascade(0), targetZ(0), targetA(0),
63 isInitialised(false)
64{
65 SetProcessSubType(fHadronAtRest);
66 Cascade = new G4GHEKinematicsVector [17];
67 pSelector = new G4StopElementSelector();
68 pEMCascade = new G4MuMinusCaptureCascade();
69 theN = new G4Fancy3DNucleus();
70 theHandler = new G4ExcitationHandler();
71 G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
72}
73
74//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75
76G4MuonMinusCaptureAtRest::~G4MuonMinusCaptureAtRest()
77{
78 G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this);
79 delete [] Cascade;
80 delete pSelector;
81 delete pEMCascade;
82 delete theN;
83 delete theHandler;
84}
85
86//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
87
88G4bool G4MuonMinusCaptureAtRest::IsApplicable(const G4ParticleDefinition& p)
89{
90 return ( &p == G4MuonMinus::MuonMinus() );
91}
92
93//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
94
95void G4MuonMinusCaptureAtRest::PreparePhysicsTable(const G4ParticleDefinition& p)
96{
97 G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this, &p);
98}
99
100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101
102void G4MuonMinusCaptureAtRest::BuildPhysicsTable(const G4ParticleDefinition& p)
103{
104 G4HadronicProcessStore::Instance()->PrintInfo(&p);
105}
106
107//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
108
109G4VParticleChange* 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
223G4ReactionProductVector* 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
Note: See TracBrowser for help on using the repository browser.