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

Last change on this file since 846 was 819, checked in by garnier, 16 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.