G4MuMinusCaptureCascade.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 //   G4MuonMinusCaptureAtRest physics process
00029 //                   
00030 //   E-mail: Vladimir.Ivantchenko@cern.ch
00031 //
00032 //   Created:   02.04.00 V.Ivanchenko
00033 //
00034 // Modified:  
00035 // 06.04.01 V.Ivanchenko Bug in theta distribution fixed
00036 // 13.02.07 V.Ivanchenko Fixes in decay - add random distribution of e- 
00037 //                       direction; factor 2 in potential energy 
00038 //
00039 //----------------------------------------------------------------------
00040 
00041 #include "G4MuMinusCaptureCascade.hh"
00042 #include "G4PhysicalConstants.hh"
00043 #include "G4SystemOfUnits.hh"
00044 #include "G4LorentzVector.hh"
00045 #include "G4ParticleMomentum.hh"
00046 #include "G4MuonMinus.hh"
00047 #include "G4Electron.hh"
00048 #include "G4Gamma.hh"
00049 #include "G4NeutrinoMu.hh"
00050 #include "G4AntiNeutrinoE.hh"
00051 #include "G4GHEKinematicsVector.hh"
00052 
00053 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00054 
00055 G4MuMinusCaptureCascade::G4MuMinusCaptureCascade()
00056 { 
00057   theElectron = G4Electron::Electron();
00058   theGamma = G4Gamma::Gamma();
00059   Emass = theElectron->GetPDGMass();
00060   MuMass = G4MuonMinus::MuonMinus()->GetPDGMass();
00061 }
00062 
00063 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00064 
00065 G4MuMinusCaptureCascade::~G4MuMinusCaptureCascade()
00066 { }
00067 
00068 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00069 
00070 G4double G4MuMinusCaptureCascade::GetKShellEnergy(G4double Z)
00071 { 
00072   // Calculate the Energy of K Mesoatom Level for this Element using
00073   // the Energy of Hydrogen Atom taken into account finite size of the
00074   // nucleus (V.Ivanchenko)
00075   const G4int ListK = 28;
00076   static G4double ListZK[ListK] = {
00077       1., 2.,  4.,  6.,  8., 11., 14., 17., 18., 21., 24.,
00078      26., 29., 32., 38., 40., 41., 44., 49., 53., 55.,
00079      60., 65., 70., 75., 81., 85., 92.};
00080   static G4double ListKEnergy[ListK] = {
00081      0.00275, 0.011, 0.043, 0.098, 0.173, 0.326,
00082      0.524, 0.765, 0.853, 1.146, 1.472,
00083      1.708, 2.081, 2.475, 3.323, 3.627, 
00084      3.779, 4.237, 5.016, 5.647, 5.966,
00085      6.793, 7.602, 8.421, 9.249, 10.222,
00086     10.923,11.984};
00087 
00088   // Energy with finit size corrections
00089   G4double KEnergy = GetLinApprox(ListK,ListZK,ListKEnergy,Z);
00090 
00091   return KEnergy;
00092 }
00093 
00094 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00095 
00096 void G4MuMinusCaptureCascade::AddNewParticle(G4ParticleDefinition* aParticle,
00097                                              G4ThreeVector& Momentum,
00098                                              G4double mass,
00099                                              G4int* nParticle,
00100                                              G4GHEKinematicsVector* Cascade)
00101 {
00102   // Store particle in the HEK vector and increment counter
00103   Cascade[*nParticle].SetZero();
00104   Cascade[*nParticle].SetMass( mass );
00105   Cascade[*nParticle].SetMomentumAndUpdate(Momentum.x(), Momentum.y(), Momentum.z());
00106   Cascade[*nParticle].SetParticleDef( aParticle );
00107   (*nParticle)++;
00108 
00109   return;
00110 }
00111 
00112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00113 
00114 G4int G4MuMinusCaptureCascade::DoCascade(const G4double Z, const G4double massA, 
00115                                          G4GHEKinematicsVector* Cascade)
00116 {
00117   // Inicialization - cascade start from 14th level
00118   // N.C.Mukhopadhyay Phy. Rep. 30 (1977) 1.
00119   G4int nPart = 0;
00120   G4double EnergyLevel[14];
00121 
00122   G4double mass = MuMass * massA / (MuMass + massA) ;
00123 
00124   const G4double KEnergy = 13.6 * eV * Z * Z * mass/ electron_mass_c2;
00125 
00126   EnergyLevel[0] = GetKShellEnergy(Z);
00127   for( G4int i = 2; i < 15; i++ ) {
00128     EnergyLevel[i-1] = KEnergy / (i*i) ;
00129   }
00130 
00131   G4int nElec  = G4int(Z);
00132   G4int nAuger = 1;
00133   G4int nLevel = 13;
00134   G4double DeltaE;
00135   G4double pGamma = Z*Z*Z*Z;
00136 
00137   // Capture on 14-th level
00138   G4double ptot = std::sqrt(EnergyLevel[13]*(EnergyLevel[13] + 2.0*Emass));
00139   G4ThreeVector moment = ptot * GetRandomVec();
00140 
00141   AddNewParticle(theElectron,moment,Emass,&nPart,Cascade);
00142 
00143   // Emit new photon or electron
00144   // Simplified model for probabilities
00145   // N.C.Mukhopadhyay Phy. Rep. 30 (1977) 1.
00146   do {
00147 
00148     // case of Auger electrons
00149     if((nAuger < nElec) && ((pGamma + 10000.0) * G4UniformRand() < 10000.0) ) {
00150         nAuger++;
00151         DeltaE = EnergyLevel[nLevel-1] - EnergyLevel[nLevel];
00152         nLevel--;
00153 
00154         ptot = std::sqrt(DeltaE * (DeltaE + 2.0*Emass));
00155         moment = ptot * GetRandomVec();
00156 
00157         AddNewParticle(theElectron, moment, Emass, &nPart, Cascade);
00158 
00159     } else {
00160 
00161       // Case of photon cascade, probabilities from
00162       // C.S.Wu and L.Wilets, Ann. Rev. Nuclear Sci. 19 (1969) 527.
00163 
00164       G4double var = (10.0 + G4double(nLevel - 1) ) * G4UniformRand();
00165       G4int iLevel = nLevel - 1 ;
00166       if(var > 10.0) iLevel -= G4int(var-10.0) + 1;
00167       if( iLevel < 0 ) iLevel = 0;
00168       DeltaE = EnergyLevel[iLevel] - EnergyLevel[nLevel];
00169       nLevel = iLevel;
00170       moment = DeltaE * GetRandomVec();
00171       AddNewParticle(theGamma, moment, 0.0, &nPart, Cascade);
00172     }
00173 
00174   } while( nLevel > 0 );
00175 
00176   return nPart;
00177 }
00178 
00179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00180 
00181 void G4MuMinusCaptureCascade::DoBoundMuonMinusDecay(G4double Z, 
00182                                                     G4int* nCascade, 
00183                                                     G4GHEKinematicsVector* Cascade)
00184 {
00185   // Simulation on Decay of mu- on a K-shell of the muonic atom
00186   G4double xmax = ( 1.0 + Emass*Emass/ (MuMass*MuMass) );
00187   G4double xmin = 2.0*Emass/MuMass;
00188   G4double KEnergy = GetKShellEnergy(Z);
00189   /*
00190   G4cout << "G4MuMinusCaptureCascade::DoBoundMuonMinusDecay" 
00191          << " XMAX= " << xmax
00192          << " Ebound= " << KEnergy
00193          << G4endl;
00194   */
00195   G4double pmu = std::sqrt(KEnergy*(KEnergy + 2.0*MuMass));
00196   G4double emu = KEnergy + MuMass;
00197   G4ThreeVector moment = GetRandomVec();
00198   G4LorentzVector MU(pmu*moment,emu);
00199   G4ThreeVector bst = MU.boostVector();
00200 
00201   G4double Eelect, Pelect, x, ecm;
00202   G4LorentzVector EL, NN;
00203   // Calculate electron energy
00204   do {
00205     do {
00206       x = xmin + (xmax-xmin)*G4UniformRand();
00207     } while (G4UniformRand() > (3.0 - 2.0*x)*x*x );
00208     Eelect = x*MuMass*0.5;
00209     Pelect = 0.0;
00210     if(Eelect > Emass) { 
00211       Pelect = std::sqrt( Eelect*Eelect - Emass*Emass );
00212     } else {
00213       Pelect = 0.0;
00214       Eelect = Emass;
00215     }
00216     G4ThreeVector e_mom = GetRandomVec();
00217     EL = G4LorentzVector(Pelect*e_mom,Eelect);
00218     EL.boost(bst);
00219     Eelect = EL.e() - Emass - 2.0*KEnergy;
00220     //
00221     // Calculate rest frame parameters of 2 neutrinos
00222     //
00223     NN = MU - EL;
00224     ecm = NN.mag2();
00225   } while (Eelect < 0.0 || ecm < 0.0);
00226 
00227   //
00228   // Create electron
00229   //
00230   moment = std::sqrt(Eelect * (Eelect + 2.0*Emass))*(EL.vect().unit());
00231   AddNewParticle(theElectron, moment, Emass, nCascade, Cascade);
00232   //
00233   // Create Neutrinos
00234   //
00235   ecm = 0.5*std::sqrt(ecm);
00236   bst = NN.boostVector();
00237   G4ThreeVector p1 = ecm * GetRandomVec();
00238   G4LorentzVector N1 = G4LorentzVector(p1,ecm);
00239   N1.boost(bst);
00240   G4ThreeVector p1lab = N1.vect();
00241   AddNewParticle(G4AntiNeutrinoE::AntiNeutrinoE(),p1lab,0.0,nCascade,Cascade);
00242   NN -= N1;
00243   G4ThreeVector p2lab = NN.vect();
00244   AddNewParticle(G4NeutrinoMu::NeutrinoMu(),p2lab,0.0,nCascade,Cascade);
00245 
00246   return;
00247 }
00248 
00249 
00250 
00251 
00252 
00253 
00254 
00255 
00256 
00257 
00258 
00259 
00260 
00261 
00262 
00263 

Generated on Mon May 27 17:48:54 2013 for Geant4 by  doxygen 1.4.7