G4PenelopeAnnihilationModel.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 // Author: Luciano Pandola
00029 //
00030 // History:
00031 // --------
00032 // 29 Oct 2008   L Pandola    Migration from process to model 
00033 // 15 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
00034 //                  - apply internal high-energy limit only in constructor 
00035 //                  - do not apply low-energy limit (default is 0)
00036 //                  - do not use G4ElementSelector
00037 
00038 #include "G4PenelopeAnnihilationModel.hh"
00039 #include "G4PhysicalConstants.hh"
00040 #include "G4SystemOfUnits.hh"
00041 #include "G4ParticleDefinition.hh"
00042 #include "G4MaterialCutsCouple.hh"
00043 #include "G4ProductionCutsTable.hh"
00044 #include "G4DynamicParticle.hh"
00045 #include "G4Gamma.hh"
00046 
00047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00048 
00049 
00050 G4PenelopeAnnihilationModel::G4PenelopeAnnihilationModel(const G4ParticleDefinition*,
00051                                              const G4String& nam)
00052   :G4VEmModel(nam),fParticleChange(0),isInitialised(false)
00053 {
00054   fIntrinsicLowEnergyLimit = 0.0;
00055   fIntrinsicHighEnergyLimit = 100.0*GeV;
00056   //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
00057   SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
00058  
00059   //Calculate variable that will be used later on
00060   fPielr2 = pi*classic_electr_radius*classic_electr_radius;
00061 
00062   verboseLevel= 0;
00063   // Verbosity scale:
00064   // 0 = nothing 
00065   // 1 = warning for energy non-conservation 
00066   // 2 = details of energy budget
00067   // 3 = calculation of cross sections, file openings, sampling of atoms
00068   // 4 = entering in methods
00069 
00070 }
00071 
00072 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00073 
00074 G4PenelopeAnnihilationModel::~G4PenelopeAnnihilationModel()
00075 {;}
00076 
00077 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00078 
00079 void G4PenelopeAnnihilationModel::Initialise(const G4ParticleDefinition*,
00080                                              const G4DataVector&)
00081 {
00082   if (verboseLevel > 3)
00083     G4cout << "Calling G4PenelopeAnnihilationModel::Initialise()" << G4endl;
00084 
00085   if(verboseLevel > 0) {
00086     G4cout << "Penelope Annihilation model is initialized " << G4endl
00087            << "Energy range: "
00088            << LowEnergyLimit() / keV << " keV - "
00089            << HighEnergyLimit() / GeV << " GeV"
00090            << G4endl;
00091   }
00092 
00093   if(isInitialised) return;
00094   fParticleChange = GetParticleChangeForGamma();
00095   isInitialised = true; 
00096 }
00097 
00098 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00099 
00100 G4double G4PenelopeAnnihilationModel::ComputeCrossSectionPerAtom(
00101                                        const G4ParticleDefinition*,
00102                                              G4double energy,
00103                                              G4double Z, G4double,
00104                                              G4double, G4double)
00105 {
00106   if (verboseLevel > 3)
00107     G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopeAnnihilationModel" << 
00108       G4endl;
00109 
00110   G4double cs = Z*ComputeCrossSectionPerElectron(energy);
00111   
00112   if (verboseLevel > 2)
00113     G4cout << "Annihilation cross Section at " << energy/keV << " keV for Z=" << Z << 
00114       " = " << cs/barn << " barn" << G4endl;
00115   return cs;
00116 }
00117 
00118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00119 
00120 void G4PenelopeAnnihilationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00121                                               const G4MaterialCutsCouple*,
00122                                               const G4DynamicParticle* aDynamicPositron,
00123                                               G4double,
00124                                               G4double)
00125 {
00126   //
00127   // Penelope model to sample final state for positron annihilation. 
00128   // Target eletrons are assumed to be free and at rest. Binding effects enabling 
00129   // one-photon annihilation are neglected.
00130   // For annihilation at rest, two back-to-back photons are emitted, having energy of 511 keV 
00131   // and isotropic angular distribution.
00132   // For annihilation in flight, it is used the theory from 
00133   //  W. Heitler, The quantum theory of radiation, Oxford University Press (1954)
00134   // The two photons can have different energy. The efficiency of the sampling algorithm 
00135   // of the photon energy from the dSigma/dE distribution is practically 100% for 
00136   // positrons of kinetic energy < 10 keV. It reaches a minimum (about 80%) at energy 
00137   // of about 10 MeV.
00138   // The angle theta is kinematically linked to the photon energy, to ensure momentum 
00139   // conservation. The angle phi is sampled isotropically for the first gamma.
00140   //
00141   if (verboseLevel > 3)
00142     G4cout << "Calling SamplingSecondaries() of G4PenelopeAnnihilationModel" << G4endl;
00143 
00144   G4double kineticEnergy = aDynamicPositron->GetKineticEnergy();
00145 
00146   // kill primary
00147   fParticleChange->SetProposedKineticEnergy(0.);
00148   fParticleChange->ProposeTrackStatus(fStopAndKill);
00149   
00150   if (kineticEnergy == 0.0)
00151     {
00152       //Old AtRestDoIt
00153       G4double cosTheta = -1.0+2.0*G4UniformRand();
00154       G4double sinTheta = std::sqrt(1.0-cosTheta*cosTheta);
00155       G4double phi = twopi*G4UniformRand();
00156       G4ThreeVector direction (sinTheta*std::cos(phi),sinTheta*std::sin(phi),cosTheta);
00157       G4DynamicParticle* firstGamma = new G4DynamicParticle (G4Gamma::Gamma(),
00158                                                              direction, electron_mass_c2);
00159       G4DynamicParticle* secondGamma = new G4DynamicParticle (G4Gamma::Gamma(),
00160                                                               -direction, electron_mass_c2);
00161   
00162       fvect->push_back(firstGamma);
00163       fvect->push_back(secondGamma);
00164       return;
00165     }
00166 
00167   //This is the "PostStep" case (annihilation in flight)
00168   G4ParticleMomentum positronDirection = 
00169     aDynamicPositron->GetMomentumDirection();
00170   G4double gamma = 1.0 + std::max(kineticEnergy,1.0*eV)/electron_mass_c2;
00171   G4double gamma21 = std::sqrt(gamma*gamma-1);
00172   G4double ani = 1.0+gamma;
00173   G4double chimin = 1.0/(ani+gamma21);
00174   G4double rchi = (1.0-chimin)/chimin;
00175   G4double gt0 = ani*ani-2.0;
00176   G4double test=0.0;
00177   G4double epsilon = 0;
00178   do{
00179     epsilon = chimin*std::pow(rchi,G4UniformRand());
00180     G4double reject = ani*ani*(1.0-epsilon)+2.0*gamma-(1.0/epsilon);
00181     test = G4UniformRand()*gt0-reject;
00182   }while(test>0);
00183    
00184   G4double totalAvailableEnergy = kineticEnergy + 2.0*electron_mass_c2;
00185   G4double photon1Energy = epsilon*totalAvailableEnergy;
00186   G4double photon2Energy = (1.0-epsilon)*totalAvailableEnergy;
00187   G4double cosTheta1 = (ani-1.0/epsilon)/gamma21;
00188   G4double cosTheta2 = (ani-1.0/(1.0-epsilon))/gamma21;
00189   
00190   //G4double localEnergyDeposit = 0.; 
00191 
00192   G4double sinTheta1 = std::sqrt(1.-cosTheta1*cosTheta1);
00193   G4double phi1  = twopi * G4UniformRand();
00194   G4double dirx1 = sinTheta1 * std::cos(phi1);
00195   G4double diry1 = sinTheta1 * std::sin(phi1);
00196   G4double dirz1 = cosTheta1;
00197   
00198   G4double sinTheta2 = std::sqrt(1.-cosTheta2*cosTheta2);
00199   G4double phi2  = phi1+pi;
00200   G4double dirx2 = sinTheta2 * std::cos(phi2);
00201   G4double diry2 = sinTheta2 * std::sin(phi2);
00202   G4double dirz2 = cosTheta2;
00203   
00204   G4ThreeVector photon1Direction (dirx1,diry1,dirz1);
00205   photon1Direction.rotateUz(positronDirection);   
00206   // create G4DynamicParticle object for the particle1  
00207   G4DynamicParticle* aParticle1= new G4DynamicParticle (G4Gamma::Gamma(),
00208                                                            photon1Direction, 
00209                                                            photon1Energy);
00210   fvect->push_back(aParticle1);
00211  
00212   G4ThreeVector photon2Direction(dirx2,diry2,dirz2);
00213   photon2Direction.rotateUz(positronDirection); 
00214      // create G4DynamicParticle object for the particle2 
00215   G4DynamicParticle* aParticle2= new G4DynamicParticle (G4Gamma::Gamma(),
00216                                                            photon2Direction,
00217                                                            photon2Energy);
00218   fvect->push_back(aParticle2);
00219 
00220   if (verboseLevel > 1)
00221     {
00222       G4cout << "-----------------------------------------------------------" << G4endl;
00223       G4cout << "Energy balance from G4PenelopeAnnihilation" << G4endl;
00224       G4cout << "Kinetic positron energy: " << kineticEnergy/keV << " keV" << G4endl;
00225       G4cout << "Total available energy: " << totalAvailableEnergy/keV << " keV " << G4endl;
00226       G4cout << "-----------------------------------------------------------" << G4endl;
00227       G4cout << "Photon energy 1: " << photon1Energy/keV << " keV" << G4endl;
00228       G4cout << "Photon energy 2: " << photon2Energy/keV << " keV" << G4endl;
00229       G4cout << "Total final state: " << (photon1Energy+photon2Energy)/keV << 
00230         " keV" << G4endl;
00231       G4cout << "-----------------------------------------------------------" << G4endl;
00232     }
00233   if (verboseLevel > 0)
00234     {      
00235       G4double energyDiff = std::fabs(totalAvailableEnergy-photon1Energy-photon2Energy);
00236       if (energyDiff > 0.05*keV)
00237         G4cout << "Warning from G4PenelopeAnnihilation: problem with energy conservation: " << 
00238           (photon1Energy+photon2Energy)/keV << 
00239           " keV (final) vs. " << 
00240           totalAvailableEnergy/keV << " keV (initial)" << G4endl;
00241     }
00242   return;
00243 }
00244 
00245 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00246 
00247 G4double G4PenelopeAnnihilationModel:: ComputeCrossSectionPerElectron(G4double energy)
00248 {
00249   //
00250   // Penelope model to calculate cross section for positron annihilation.
00251   // The annihilation cross section per electron is calculated according 
00252   // to the Heitler formula
00253   //  W. Heitler, The quantum theory of radiation, Oxford University Press (1954)
00254   // in the assumptions of electrons free and at rest.
00255   //
00256   G4double gamma = 1.0+std::max(energy,1.0*eV)/electron_mass_c2;
00257   G4double gamma2 = gamma*gamma;
00258   G4double f2 = gamma2-1.0;
00259   G4double f1 = std::sqrt(f2);
00260   G4double crossSection = fPielr2*((gamma2+4.0*gamma+1.0)*std::log(gamma+f1)/f2
00261                          - (gamma+3.0)/f1)/(gamma+1.0);
00262   return crossSection;
00263 }

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