Geant4_10
G4MuonMinusCaptureAtRest.cc
Go to the documentation of this file.
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 70681 2013-06-04 07:54:04Z gcosmo $
27 //
28 // G4MuonMinusCaptureAtRest physics process
29 // Larry Felawka (TRIUMF) and Art Olin (TRIUMF)
30 // April 1998
31 //---------------------------------------------------------------------
32 //
33 // Modifications:
34 // 18/08/2000 V.Ivanchenko Update description
35 // 12/12/2003 H.P.Wellisch Completly rewrite mu-nuclear part
36 // 17/05/2006 V.Ivanchenko Cleanup
37 // 15/11/2006 V.Ivanchenko Review and rewrite all kinematics
38 // 24/01/2007 V.Ivanchenko Force to work with integer Z and A
39 // 23/01/2009 V.Ivanchenko Add deregistration
40 //
41 //-----------------------------------------------------------------------------
42 
44 #include "G4DynamicParticle.hh"
45 #include "Randomize.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4He3.hh"
49 #include "G4NeutrinoMu.hh"
50 #include "G4Fragment.hh"
52 #include "G4Proton.hh"
53 #include "G4PionPlus.hh"
54 #include "G4GHEKinematicsVector.hh"
55 #include "G4Fancy3DNucleus.hh"
56 #include "G4ExcitationHandler.hh"
58 
59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
60 
62  G4ProcessType aType ) :
63  G4VRestProcess (processName, aType), nCascade(0), targetZ(0), targetA(0),
64  isInitialised(false)
65 {
67  Cascade = new G4GHEKinematicsVector [17];
68  pSelector = new G4StopElementSelector();
69  pEMCascade = new G4MuMinusCaptureCascade();
70  theN = new G4Fancy3DNucleus();
71  theHandler = new G4ExcitationHandler();
73 }
74 
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76 
78 {
80  delete [] Cascade;
81  delete pSelector;
82  delete pEMCascade;
83  delete theN;
84  delete theHandler;
85 }
86 
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88 
90 {
91  return ( &p == G4MuonMinus::MuonMinus() );
92 }
93 
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
95 
97 {
99 }
100 
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102 
104 {
106 }
107 
108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
109 
111  const G4Step&)
112 {
113  //
114  // Handles MuonMinuss at rest; a MuonMinus can either create secondaries or
115  // do nothing (in which case it should be sent back to decay-handling
116  // section
117  //
119 
120  // select element and get Z,A.
121  G4Element* aEle = pSelector->GetElement(track.GetMaterial());
122  targetZ = aEle->GetZ();
123  targetA = G4double(G4int(aEle->GetN()+0.5));
124  G4int ni = 0;
125 
126  G4IsotopeVector* isv = aEle->GetIsotopeVector();
127  if(isv) ni = isv->size();
128 
129  if(ni == 1) {
130  targetA = G4double(aEle->GetIsotope(0)->GetN());
131  } else if(ni > 1) {
132  G4double* ab = aEle->GetRelativeAbundanceVector();
134  G4int j = -1;
135  ni--;
136  do {
137  j++;
138  y -= ab[j];
139  } while (y > 0.0 && j < ni);
140  targetA = G4double(aEle->GetIsotope(j)->GetN());
141  }
142 
143  // Do the electromagnetic cascade of the muon in the nuclear field.
144  nCascade = 0;
145  targetMass = G4NucleiProperties::GetNuclearMass(targetA, targetZ);
146  nCascade = pEMCascade->DoCascade(targetZ, targetMass, Cascade);
147 
148  // Decide on Decay or Capture, and doit.
149  G4double lambdac = pSelector->GetMuonCaptureRate(targetZ, targetA);
150  G4double lambdad = pSelector->GetMuonDecayRate(targetZ, targetA);
151  G4double lambda = lambdac + lambdad;
152 
153  // === Throw for capture time.
154 
155  G4double tDelay = -std::log(G4UniformRand()) / lambda;
156 
157  G4ReactionProductVector * captureResult = 0;
158  G4int nEmSecondaries = nCascade;
159  G4int nSecondaries = nCascade;
160  /*
161  G4cout << "lambda= " << lambda << " lambdac= " << lambdac
162  << " nem= " << nEmSecondaries << G4endl;
163  */
164  if( G4UniformRand()*lambda > lambdac)
165  pEMCascade->DoBoundMuonMinusDecay(targetZ, &nEmSecondaries, Cascade);
166  else
167  captureResult = DoMuCapture();
168 
169  // fill the final state
170  if(captureResult) nSecondaries += captureResult->size();
171  else nSecondaries = nEmSecondaries;
172  //G4cout << " nsec= " << nSecondaries << " nem= " << nEmSecondaries << G4endl;
173 
174  aParticleChange.SetNumberOfSecondaries( nSecondaries );
175 
176  G4double globalTime = track.GetGlobalTime();
178  // Store nuclear cascade
179  if(captureResult) {
180  G4int n = captureResult->size();
181  for ( G4int isec = 0; isec < n; isec++ ) {
182  G4ReactionProduct* aParticle = captureResult->operator[](isec);
183  G4DynamicParticle * aNewParticle = new G4DynamicParticle();
184  aNewParticle->SetDefinition( aParticle->GetDefinition() );
185  G4LorentzVector itV(aParticle->GetTotalEnergy(), aParticle->GetMomentum());
186  aNewParticle->SetMomentum(itV.vect());
187  G4double localtime = globalTime + tDelay + aParticle->GetTOF();
188  G4Track* aNewTrack = new G4Track( aNewParticle, localtime, position);
189  aNewTrack->SetTouchableHandle(track.GetTouchableHandle());
190  aParticleChange.AddSecondary( aNewTrack );
191  delete aParticle;
192  }
193  delete captureResult;
194  }
195 
196  // Store electromagnetic cascade
197 
198  if(nEmSecondaries > 0) {
199 
200  for ( G4int isec = 0; isec < nEmSecondaries; isec++ ) {
201  G4ParticleDefinition* pd = Cascade[isec].GetParticleDef();
202  G4double localtime = globalTime;
203  if(isec >= nCascade) localtime += tDelay;
204  if(pd) {
205  G4DynamicParticle* aNewParticle = new G4DynamicParticle;
206  aNewParticle->SetDefinition( pd );
207  aNewParticle->SetMomentum( Cascade[isec].GetMomentum() );
208 
209  G4Track* aNewTrack = new G4Track( aNewParticle, localtime, position );
210  aNewTrack->SetTouchableHandle(track.GetTouchableHandle());
211  aParticleChange.AddSecondary( aNewTrack );
212  }
213  }
214  }
215 
218 
219  return &aParticleChange;
220 }
221 
222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
223 
224 G4ReactionProductVector* G4MuonMinusCaptureAtRest::DoMuCapture()
225 {
227  G4double muBindingEnergy = pEMCascade->GetKShellEnergy(targetZ);
228  /*
229  G4cout << "G4MuonMinusCaptureAtRest::DoMuCapture called Emu= "
230  << muBindingEnergy << G4endl;
231  */
232  // Energy on K-shell
233  G4double muEnergy = mumass + muBindingEnergy;
234  G4double muMom = std::sqrt(muBindingEnergy*(muBindingEnergy + 2.0*mumass));
235  G4double availableEnergy = targetMass + mumass - muBindingEnergy;
236  G4LorentzVector momInitial(0.0,0.0,0.0,availableEnergy);
237  G4LorentzVector momResidual;
238 
239  G4ThreeVector vmu = muMom*pEMCascade->GetRandomVec();
240  G4LorentzVector aMuMom(vmu, muEnergy);
241 
242  G4double residualMass =
243  G4NucleiProperties::GetNuclearMass(targetA, targetZ - 1.0);
244 
245  G4ReactionProductVector* aPreResult = 0;
248 
249  G4int iz = G4int(targetZ);
250  G4int ia = G4int(targetA);
251 
252  // p, d, t, 3He or alpha as target
253  if(iz <= 2) {
254 
255  if(ia > 1) {
256  if(iz == 1 && ia == 2) {
257  availableEnergy -= neutron_mass_c2;
258  } else if(iz == 1 && ia == 3) {
259  availableEnergy -= 2.0*neutron_mass_c2;
260  } else if(iz == 2) {
261  G4ParticleDefinition* pd = 0;
262  if (ia == 3) {
263  pd = G4Deuteron::Deuteron();
264  } else if(ia == 4) {
265  pd = G4Triton::Triton();
266  } else {
268  }
269 
270  // G4cout << "Extra " << pd->GetParticleName() << G4endl;
271  availableEnergy -= pd->GetPDGMass();
272  }
273  }
274  //
275  // Computation in assumption of CM collision of mu and nucleaon
276  //
277  G4double Enu = 0.5*(availableEnergy -
278  neutron_mass_c2*neutron_mass_c2/availableEnergy);
279 
280  // make the nu, and transform to lab;
281  G4ThreeVector nu3Mom = Enu*pEMCascade->GetRandomVec();
282 
285  aN->SetTotalEnergy( availableEnergy - Enu );
286  aN->SetMomentum( -nu3Mom );
287 
288  aNu->SetTotalEnergy( Enu );
289  aNu->SetMomentum( nu3Mom );
290  aPreResult = new G4ReactionProductVector();
291 
292  aPreResult->push_back(aN );
293  aPreResult->push_back(aNu);
294 
295  if(verboseLevel > 1)
296  G4cout << "DoMuCapture on H or He"
297  <<" EkinN(MeV)= " << (availableEnergy - Enu - neutron_mass_c2)/MeV
298  <<" Enu(MeV)= "<<aNu->GetTotalEnergy()/MeV
299  <<" n= " << aPreResult->size()
300  <<G4endl;
301 
302  return aPreResult;
303  }
304 
305  // pick random proton inside nucleus
306  G4double eEx;
307  do {
308  theN->Init(ia, iz);
309  G4LorentzVector thePMom;
310  G4Nucleon * aNucleon = 0;
311  G4int theProtonCounter = G4int( targetZ * G4UniformRand() );
312  G4int counter = 0;
313  theN->StartLoop();
314 
315  while( (aNucleon=theN->GetNextNucleon()) ) {
316 
317  if( aNucleon->GetDefinition() == G4Proton::Proton() ) {
318  counter++;
319  if(counter == theProtonCounter) {
320  thePMom = aNucleon->GetMomentum();
321  break;
322  }
323  }
324  }
325 
326  // Get the nu momentum in the CMS
327  G4LorentzVector theCMS = thePMom + aMuMom;
328  G4ThreeVector bst = theCMS.boostVector();
329 
330  G4double Ecms = theCMS.mag();
331  G4double Enu = 0.5*(Ecms - neutron_mass_c2*neutron_mass_c2/Ecms);
332  eEx = 0.0;
333 
334  if(Enu > 0.0) {
335  // make the nu, and transform to lab;
336  G4ThreeVector nu3Mom = Enu*pEMCascade->GetRandomVec();
337  G4LorentzVector nuMom(nu3Mom, Enu);
338 
339  // nu in lab.
340  nuMom.boost(bst);
341  aNu->SetTotalEnergy( nuMom.e() );
342  aNu->SetMomentum( nuMom.vect() );
343 
344  // make residual
345  momResidual = momInitial - nuMom;
346 
347  // Call pre-compound on the rest.
348  eEx = momResidual.mag();
349  if(verboseLevel > 1)
350  G4cout << "G4MuonMinusCaptureAtRest::DoMuCapture: "
351  << " Eex(MeV)= " << (eEx-residualMass)/MeV
352  << " Enu(MeV)= "<<aNu->GetTotalEnergy()/MeV
353  <<G4endl;
354  }
355  } while(eEx <= residualMass);
356 
357 // G4cout << "muonCapture : " << eEx << " " << residualMass
358 // << " A,Z= " << targetA << ", "<< targetZ
359 // << " " << G4int(targetA) << ", " << G4int(targetZ) << G4endl;
360 
361  //
362  // Start Deexcitation
363  //
364  G4ThreeVector fromBreit = momResidual.boostVector();
365  G4LorentzVector fscm(0.0,0.0,0.0, eEx);
366  G4Fragment anInitialState;
367  anInitialState.SetA(targetA);
368  anInitialState.SetZ(G4double(iz - 1));
369  anInitialState.SetNumberOfParticles(2);
370  anInitialState.SetNumberOfCharged(0);
371  anInitialState.SetNumberOfHoles(1);
372  anInitialState.SetMomentum(fscm);
373  aPreResult = theHandler->BreakItUp(anInitialState);
374 
375  G4ReactionProductVector::iterator ires;
376  G4double eBal = availableEnergy - aNu->GetTotalEnergy();
377  for(ires=aPreResult->begin(); ires!=aPreResult->end(); ires++) {
378  G4LorentzVector itV((*ires)->GetTotalEnergy(), (*ires)->GetMomentum());
379  itV.boost(fromBreit);
380  (*ires)->SetTotalEnergy(itV.t());
381  (*ires)->SetMomentum(itV.vect());
382  eBal -= itV.t();
383  }
384  //
385  // fill neutrino into result
386  //
387  aPreResult->push_back(aNu);
388 
389  if(verboseLevel > 1)
390  G4cout << "DoMuCapture: Nsec= "
391  << aPreResult->size() << " Ebalance(MeV)= " << eBal/MeV
392  <<" E0(MeV)= " <<availableEnergy/MeV
393  <<" Mres(GeV)= " <<residualMass/GeV
394  <<G4endl;
395 
396  return aPreResult;
397 }
398 
void DeRegisterExtraProcess(G4VProcess *)
Hep3Vector boostVector() const
G4MuonMinusCaptureAtRest(const G4String &processName="muMinusCaptureAtRest", G4ProcessType aType=fHadronic)
static G4double GetNuclearMass(const G4double A, const G4double Z)
std::vector< G4Isotope * > G4IsotopeVector
void SetMomentum(const G4ThreeVector &momentum)
const G4LorentzVector & GetMomentum() const
Definition: G4Nucleon.hh:71
G4int verboseLevel
Definition: G4VProcess.hh:368
static G4HadronicProcessStore * Instance()
G4double GetN() const
Definition: G4Element.hh:134
G4ParticleDefinition * GetParticleDef()
void SetMomentum(const G4double x, const G4double y, const G4double z)
const char * p
Definition: xmltok.h:285
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:449
virtual void PreparePhysicsTable(const G4ParticleDefinition &)
const G4ThreeVector & GetPosition() const
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
G4double GetZ() const
Definition: G4Element.hh:131
G4double GetMuonCaptureRate(G4double Z, G4double A)
void SetNumberOfHoles(G4int valueTot, G4int valueP=0)
Definition: G4Fragment.hh:335
G4double GetMuonDecayRate(G4double Z, G4double A)
void SetTouchableHandle(const G4TouchableHandle &apValue)
int G4int
Definition: G4Types.hh:78
G4ReactionProductVector * BreakItUp(const G4Fragment &theInitialState) const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4ParticleDefinition * GetDefinition() const
void SetA(G4double value)
Definition: G4Fragment.hh:294
std::vector< G4ReactionProduct * > G4ReactionProductVector
Double_t y
Definition: plot.C:279
virtual G4bool IsApplicable(const G4ParticleDefinition &)
Char_t n[5]
void RegisterParticleForExtraProcess(G4VProcess *, const G4ParticleDefinition *)
G4IonTable * GetIonTable() const
G4int GetN() const
Definition: G4Isotope.hh:94
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
double mag() const
static G4NeutrinoMu * NeutrinoMu()
Definition: G4NeutrinoMu.cc:85
bool G4bool
Definition: G4Types.hh:79
G4int DoCascade(const G4double Z, const G4double A, G4GHEKinematicsVector *Cascade)
HepLorentzVector & boost(double, double, double)
G4double iz
Definition: TRTMaterials.hh:39
void SetMomentum(const G4LorentzVector &value)
Definition: G4Fragment.hh:256
virtual G4VParticleChange * AtRestDoIt(const G4Track &, const G4Step &)
void SetTotalEnergy(const G4double en)
void Init(G4int theA, G4int theZ)
static G4Triton * Triton()
Definition: G4Triton.cc:95
G4double * GetRelativeAbundanceVector() const
Definition: G4Element.hh:166
static G4Proton * Proton()
Definition: G4Proton.cc:93
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
Definition: G4Step.hh:76
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
void SetNumberOfParticles(G4int value)
Definition: G4Fragment.hh:344
G4double GetGlobalTime() const
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
const G4TouchableHandle & GetTouchableHandle() const
G4Material * GetMaterial() const
float neutron_mass_c2
Definition: hepunit.py:276
void RegisterExtraProcess(G4VProcess *)
virtual void Initialize(const G4Track &)
G4double GetTotalEnergy() const
G4IsotopeVector * GetIsotopeVector() const
Definition: G4Element.hh:162
G4double GetKShellEnergy(G4double Z)
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
void SetNumberOfSecondaries(G4int totSecondaries)
G4double GetTOF() const
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
G4Element * GetElement(const G4Material *aMaterial)
void SetDefinition(G4ParticleDefinition *aParticleDefinition)
void DoBoundMuonMinusDecay(G4double Z, G4int *nCascade, G4GHEKinematicsVector *Cascade)
void AddSecondary(G4Track *aSecondary)
G4ThreeVector GetMomentum() const
const G4Isotope * GetIsotope(G4int iso) const
Definition: G4Element.hh:169
static G4MuonMinus * MuonMinus()
Definition: G4MuonMinus.cc:100
#define G4endl
Definition: G4ios.hh:61
virtual G4ParticleDefinition * GetDefinition() const
Definition: G4Nucleon.hh:85
void SetZ(G4double value)
Definition: G4Fragment.hh:288
void SetNumberOfCharged(G4int value)
Definition: G4Fragment.hh:349
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4Nucleon * GetNextNucleon()
void PrintInfo(const G4ParticleDefinition *)
G4ProcessType