Geant4: G4AnnihiToMuPair.cc Source 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 //
00027 // $Id: G4AnnihiToMuPair.cc 66996 2013-01-29 14:50:52Z gcosmo $
00028 //
00029 //         ------------ G4AnnihiToMuPair physics process ------
00030 //         by H.Burkhardt, S. Kelner and R. Kokoulin, November 2002
00031 // -----------------------------------------------------------------------------
00032 //
00033 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......//
00034 //
00035 // 27.01.03 : first implementation (hbu)
00036 // 04.02.03 : cosmetic simplifications (mma)
00037 // 25.10.04 : migrade to new interfaces of ParticleChange (vi)
00038 //
00039 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00040 
00041 #include "G4AnnihiToMuPair.hh"
00042 
00043 #include "G4ios.hh"
00044 #include "Randomize.hh"
00045 #include "G4PhysicalConstants.hh"
00046 #include "G4SystemOfUnits.hh"
00047 
00048 #include "G4Positron.hh"
00049 #include "G4MuonPlus.hh"
00050 #include "G4MuonMinus.hh"
00051 #include "G4Material.hh"
00052 #include "G4Step.hh"
00053 
00054 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00055 
00056 using namespace std;
00057 
00058 G4AnnihiToMuPair::G4AnnihiToMuPair(const G4String& processName,
00059     G4ProcessType type):G4VDiscreteProcess (processName, type)
00060 {
00061  //e+ Energy threshold
00062  const G4double Mu_massc2 = G4MuonPlus::MuonPlus()->GetPDGMass();
00063  LowestEnergyLimit  = 2*Mu_massc2*Mu_massc2/electron_mass_c2 - electron_mass_c2;
00064  
00065  //modele ok up to 1000 TeV due to neglected Z-interference
00066  HighestEnergyLimit = 1000*TeV;
00067  
00068  CurrentSigma = 0.0;
00069  CrossSecFactor = 1.;
00070  SetProcessSubType(6);
00071 
00072 }
00073 
00074 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00075 
00076 G4AnnihiToMuPair::~G4AnnihiToMuPair() // (empty) destructor
00077 { }
00078 
00079 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00080 
00081 G4bool G4AnnihiToMuPair::IsApplicable(const G4ParticleDefinition& particle)
00082 {
00083   return ( &particle == G4Positron::Positron() );
00084 }
00085 
00086 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00087 
00088 void G4AnnihiToMuPair::BuildPhysicsTable(const G4ParticleDefinition&)
00089 // Build cross section and mean free path tables
00090 //here no tables, just calling PrintInfoDefinition
00091 {
00092   CurrentSigma = 0.0;
00093   PrintInfoDefinition();
00094 }
00095 
00096 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00097 
00098 void G4AnnihiToMuPair::SetCrossSecFactor(G4double fac)
00099 // Set the factor to artificially increase the cross section
00100 { 
00101   CrossSecFactor = fac;
00102   G4cout << "The cross section for AnnihiToMuPair is artificially "
00103          << "increased by the CrossSecFactor=" << CrossSecFactor << G4endl;
00104 }
00105 
00106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00107 
00108 G4double G4AnnihiToMuPair::ComputeCrossSectionPerAtom(G4double Epos, G4double Z)
00109 // Calculates the microscopic cross section in GEANT4 internal units.
00110 // It gives a good description from threshold to 1000 GeV
00111 {
00112   static const G4double Mmuon = G4MuonPlus::MuonPlus()->GetPDGMass();
00113   static const G4double Rmuon = elm_coupling/Mmuon; //classical particle radius
00114   static const G4double Sig0  = pi*Rmuon*Rmuon/3.;  //constant in crossSection
00115 
00116   G4double CrossSection = 0.;
00117   if (Epos < LowestEnergyLimit) return CrossSection;
00118    
00119   G4double xi = LowestEnergyLimit/Epos;
00120   G4double SigmaEl = Sig0*xi*(1.+xi/2.)*sqrt(1.-xi); // per electron
00121   CrossSection = SigmaEl*Z;         // number of electrons per atom
00122   return CrossSection;
00123 }
00124 
00125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00126 
00127 G4double G4AnnihiToMuPair::CrossSectionPerVolume(G4double PositronEnergy, 
00128                                                  const G4Material* aMaterial)
00129 {
00130   const G4ElementVector* theElementVector = aMaterial->GetElementVector();
00131   const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
00132 
00133   G4double SIGMA = 0.0;
00134 
00135   for ( size_t i=0 ; i < aMaterial->GetNumberOfElements() ; ++i )
00136   {
00137     G4double AtomicZ = (*theElementVector)[i]->GetZ();
00138     SIGMA += NbOfAtomsPerVolume[i] *
00139       ComputeCrossSectionPerAtom(PositronEnergy,AtomicZ);
00140   }
00141   return SIGMA;
00142 }
00143 
00144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00145 
00146 G4double G4AnnihiToMuPair::GetMeanFreePath(const G4Track& aTrack,
00147                                            G4double, G4ForceCondition*)
00148 
00149 // returns the positron mean free path in GEANT4 internal units
00150 
00151 {
00152   const G4DynamicParticle* aDynamicPositron = aTrack.GetDynamicParticle();
00153   G4double PositronEnergy = aDynamicPositron->GetKineticEnergy()
00154                                               +electron_mass_c2;
00155   G4Material* aMaterial = aTrack.GetMaterial();
00156   CurrentSigma = CrossSectionPerVolume(PositronEnergy, aMaterial);
00157 
00158   // increase the CrossSection by CrossSecFactor (default 1)
00159   G4double mfp = DBL_MAX;
00160   if(CurrentSigma > DBL_MIN) mfp = 1.0/(CurrentSigma*CrossSecFactor);
00161 
00162   return mfp;
00163 }
00164 
00165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00166 
00167 G4VParticleChange* G4AnnihiToMuPair::PostStepDoIt(const G4Track& aTrack,
00168                                                   const G4Step&  aStep)
00169 //
00170 // generation of e+e- -> mu+mu-
00171 //
00172 {
00173 
00174   aParticleChange.Initialize(aTrack);
00175   static const G4double Mele=electron_mass_c2;
00176   static const G4double Mmuon=G4MuonPlus::MuonPlus()->GetPDGMass();
00177 
00178   // current Positron energy and direction, return if energy too low
00179   const G4DynamicParticle *aDynamicPositron = aTrack.GetDynamicParticle();
00180   G4double Epos = aDynamicPositron->GetKineticEnergy() + Mele; 
00181 
00182   // test of cross section
00183   if(CurrentSigma*G4UniformRand() > 
00184      CrossSectionPerVolume(Epos, aTrack.GetMaterial())) 
00185     {
00186       return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
00187     }
00188 
00189   if (Epos < LowestEnergyLimit) {
00190      return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
00191   }
00192 
00193   G4ParticleMomentum PositronDirection = 
00194                                        aDynamicPositron->GetMomentumDirection();
00195   G4double xi = LowestEnergyLimit/Epos; // xi is always less than 1,
00196                                         // goes to 0 at high Epos
00197 
00198   // generate cost
00199   //
00200   G4double cost;
00201   do cost = 2.*G4UniformRand()-1.;
00202   while (2.*G4UniformRand() > 1.+xi+cost*cost*(1.-xi) ); 
00203                                                        //1+cost**2 at high Epos
00204   G4double sint = sqrt(1.-cost*cost);
00205 
00206   // generate phi
00207   //
00208   G4double phi=2.*pi*G4UniformRand();
00209 
00210   G4double Ecm   = sqrt(0.5*Mele*(Epos+Mele));
00211   G4double Pcm   = sqrt(Ecm*Ecm-Mmuon*Mmuon);
00212   G4double beta  = sqrt((Epos-Mele)/(Epos+Mele));
00213   G4double gamma = Ecm/Mele;                    // =sqrt((Epos+Mele)/(2.*Mele));
00214   G4double Pt    = Pcm*sint;
00215   
00216   // energy and momentum of the muons in the Lab
00217   //
00218   G4double EmuPlus   = gamma*(     Ecm+cost*beta*Pcm);
00219   G4double EmuMinus  = gamma*(     Ecm-cost*beta*Pcm);
00220   G4double PmuPlusZ  = gamma*(beta*Ecm+cost*     Pcm);
00221   G4double PmuMinusZ = gamma*(beta*Ecm-cost*     Pcm);
00222   G4double PmuPlusX  = Pt*cos(phi);
00223   G4double PmuPlusY  = Pt*sin(phi);
00224   G4double PmuMinusX =-Pt*cos(phi);
00225   G4double PmuMinusY =-Pt*sin(phi);
00226   // absolute momenta
00227   G4double PmuPlus  = sqrt(Pt*Pt+PmuPlusZ *PmuPlusZ );
00228   G4double PmuMinus = sqrt(Pt*Pt+PmuMinusZ*PmuMinusZ);
00229 
00230   // mu+ mu- directions for Positron in z-direction
00231   //
00232   G4ThreeVector
00233     MuPlusDirection ( PmuPlusX/PmuPlus, PmuPlusY/PmuPlus,  PmuPlusZ/PmuPlus  );
00234   G4ThreeVector
00235     MuMinusDirection(PmuMinusX/PmuMinus,PmuMinusY/PmuMinus,PmuMinusZ/PmuMinus);
00236 
00237   // rotate to actual Positron direction
00238   //
00239   MuPlusDirection.rotateUz(PositronDirection);
00240   MuMinusDirection.rotateUz(PositronDirection);
00241 
00242   aParticleChange.SetNumberOfSecondaries(2);
00243   // create G4DynamicParticle object for the particle1
00244   G4DynamicParticle* aParticle1= new G4DynamicParticle(
00245                          G4MuonPlus::MuonPlus(),MuPlusDirection,EmuPlus-Mmuon);
00246   aParticleChange.AddSecondary(aParticle1);
00247   // create G4DynamicParticle object for the particle2
00248   G4DynamicParticle* aParticle2= new G4DynamicParticle(
00249                      G4MuonMinus::MuonMinus(),MuMinusDirection,EmuMinus-Mmuon);
00250   aParticleChange.AddSecondary(aParticle2);
00251 
00252   // Kill the incident positron 
00253   //
00254   aParticleChange.ProposeEnergy(0.); 
00255   aParticleChange.ProposeTrackStatus(fStopAndKill);
00256 
00257   return &aParticleChange;
00258 }
00259 
00260 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00261 
00262 void G4AnnihiToMuPair::PrintInfoDefinition()
00263 {
00264   G4String comments ="e+e->mu+mu- annihilation, atomic e- at rest, SubType=.";
00265   G4cout << G4endl << GetProcessName() << ":  " << comments 
00266          << GetProcessSubType() << G4endl;
00267   G4cout << "        threshold at " << LowestEnergyLimit/GeV << " GeV"
00268          << " good description up to "
00269          << HighestEnergyLimit/TeV << " TeV for all Z." << G4endl;
00270 }
00271 
00272 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......