source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeAnnihilationModel.cc @ 989

Last change on this file since 989 was 968, checked in by garnier, 15 years ago

fichier ajoutes
File size: 11.3 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: G4PenelopeAnnihilationModel.cc,v 1.2 2008/12/04 14:09:36 pandola Exp $
27// GEANT4 tag $Name: geant4-09-02-ref-02 $
28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
33// 29 Oct 2008   L Pandola    Migration from process to model
34//
35
36#include "G4PenelopeAnnihilationModel.hh"
37#include "G4ParticleDefinition.hh"
38#include "G4MaterialCutsCouple.hh"
39#include "G4ProductionCutsTable.hh"
40#include "G4DynamicParticle.hh"
41#include "G4Gamma.hh"
42
43//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44
45
46G4PenelopeAnnihilationModel::G4PenelopeAnnihilationModel(const G4ParticleDefinition*,
47                                             const G4String& nam)
48  :G4VEmModel(nam),isInitialised(false)
49{
50  fIntrinsicLowEnergyLimit = 0.0*eV;
51  fIntrinsicHighEnergyLimit = 100.0*GeV;
52  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
53  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
54 
55  //Calculate variable that will be used later on
56  fPielr2 = pi*classic_electr_radius*classic_electr_radius;
57
58  verboseLevel= 0;
59  // Verbosity scale:
60  // 0 = nothing
61  // 1 = warning for energy non-conservation
62  // 2 = details of energy budget
63  // 3 = calculation of cross sections, file openings, sampling of atoms
64  // 4 = entering in methods
65
66}
67
68//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69
70G4PenelopeAnnihilationModel::~G4PenelopeAnnihilationModel()
71{;}
72
73//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74
75void G4PenelopeAnnihilationModel::Initialise(const G4ParticleDefinition* particle,
76                                       const G4DataVector& cuts)
77{
78  if (verboseLevel > 3)
79    G4cout << "Calling G4PenelopeAnnihilationModel::Initialise()" << G4endl;
80
81  InitialiseElementSelectors(particle,cuts);
82  if (LowEnergyLimit() < fIntrinsicLowEnergyLimit)
83    {
84      G4cout << "G4PenelopeAnnihilationModel: low energy limit increased from " << 
85        LowEnergyLimit()/eV << " eV to " << fIntrinsicLowEnergyLimit/eV << " eV" << G4endl;
86      SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
87    }
88  if (HighEnergyLimit() > fIntrinsicHighEnergyLimit)
89    {
90      G4cout << "G4PenelopeAnnihilationModel: high energy limit decreased from " << 
91        HighEnergyLimit()/GeV << " GeV to " << fIntrinsicHighEnergyLimit/GeV << " GeV" << G4endl;
92      SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
93    }
94
95  G4cout << "Penelope Annihilation model is initialized " << G4endl
96         << "Energy range: "
97         << LowEnergyLimit() / keV << " keV - "
98         << HighEnergyLimit() / GeV << " GeV"
99         << G4endl;
100
101  if(isInitialised) return;
102
103  if(pParticleChange)
104    fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
105  else
106    fParticleChange = new G4ParticleChangeForGamma();
107  isInitialised = true; 
108}
109
110//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111
112G4double G4PenelopeAnnihilationModel::ComputeCrossSectionPerAtom(
113                                       const G4ParticleDefinition*,
114                                             G4double energy,
115                                             G4double Z, G4double,
116                                             G4double, G4double)
117{
118  if (verboseLevel > 3)
119    G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopeAnnihilationModel" << 
120      G4endl;
121
122  G4double cs = Z*ComputeCrossSectionPerElectron(energy);
123 
124  if (verboseLevel > 2)
125    G4cout << "Annihilation cross Section at " << energy/keV << " keV for Z=" << Z << 
126      " = " << cs/barn << " barn" << G4endl;
127  return cs;
128}
129
130//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
131
132void G4PenelopeAnnihilationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
133                                              const G4MaterialCutsCouple*,
134                                              const G4DynamicParticle* aDynamicPositron,
135                                              G4double,
136                                              G4double)
137{
138  //
139  // Penelope model to sample final state for positron annihilation.
140  // Target eletrons are assumed to be free and at rest. Binding effects enabling
141  // one-photon annihilation are neglected.
142  // For annihilation at rest, two back-to-back photons are emitted, having energy of 511 keV
143  // and isotropic angular distribution.
144  // For annihilation in flight, it is used the theory from
145  //  W. Heitler, The quantum theory of radiation, Oxford University Press (1954)
146  // The two photons can have different energy. The efficiency of the sampling algorithm
147  // of the photon energy from the dSigma/dE distribution is practically 100% for
148  // positrons of kinetic energy < 10 keV. It reaches a minimum (about 80%) at energy
149  // of about 10 MeV.
150  // The angle theta is kinematically linked to the photon energy, to ensure momentum
151  // conservation. The angle phi is sampled isotropically for the first gamma.
152  //
153  if (verboseLevel > 3)
154    G4cout << "Calling SamplingSecondaries() of G4PenelopeAnnihilationModel" << G4endl;
155
156  G4double kineticEnergy = aDynamicPositron->GetKineticEnergy();
157 
158  if (kineticEnergy == 0)
159    {
160      //Old AtRestDoIt
161      G4double cosTheta = -1.0+2.0*G4UniformRand();
162      G4double sinTheta = std::sqrt(1.0-cosTheta*cosTheta);
163      G4double phi = twopi*G4UniformRand();
164      G4ThreeVector direction (sinTheta*std::cos(phi),sinTheta*std::sin(phi),cosTheta);
165      G4DynamicParticle* firstGamma = new G4DynamicParticle (G4Gamma::Gamma(),
166                                                             direction, electron_mass_c2);
167      G4DynamicParticle* secondGamma = new G4DynamicParticle (G4Gamma::Gamma(),
168                                                              -direction, electron_mass_c2);
169 
170      fvect->push_back(firstGamma);
171      fvect->push_back(secondGamma);
172      fParticleChange->SetProposedKineticEnergy(0.);
173      fParticleChange->ProposeTrackStatus(fStopAndKill);
174      return;
175    }
176
177  //This is the "PostStep" case (annihilation in flight)
178  G4ParticleMomentum positronDirection = 
179    aDynamicPositron->GetMomentumDirection();
180  G4double gamma = 1.0 + std::max(kineticEnergy,1.0*eV)/electron_mass_c2;
181  G4double gamma21 = std::sqrt(gamma*gamma-1);
182  G4double ani = 1.0+gamma;
183  G4double chimin = 1.0/(ani+gamma21);
184  G4double rchi = (1.0-chimin)/chimin;
185  G4double gt0 = ani*ani-2.0;
186  G4double test=0.0;
187  G4double epsilon = 0;
188  do{
189    epsilon = chimin*std::pow(rchi,G4UniformRand());
190    G4double reject = ani*ani*(1.0-epsilon)+2.0*gamma-(1.0/epsilon);
191    test = G4UniformRand()*gt0-reject;
192  }while(test>0);
193   
194  G4double totalAvailableEnergy = kineticEnergy + 2.0*electron_mass_c2;
195  G4double photon1Energy = epsilon*totalAvailableEnergy;
196  G4double photon2Energy = (1.0-epsilon)*totalAvailableEnergy;
197  G4double cosTheta1 = (ani-1.0/epsilon)/gamma21;
198  G4double cosTheta2 = (ani-1.0/(1.0-epsilon))/gamma21;
199 
200  //G4double localEnergyDeposit = 0.;
201
202  G4double sinTheta1 = std::sqrt(1.-cosTheta1*cosTheta1);
203  G4double phi1  = twopi * G4UniformRand();
204  G4double dirx1 = sinTheta1 * std::cos(phi1);
205  G4double diry1 = sinTheta1 * std::sin(phi1);
206  G4double dirz1 = cosTheta1;
207 
208  G4double sinTheta2 = std::sqrt(1.-cosTheta2*cosTheta2);
209  G4double phi2  = phi1+pi;
210  G4double dirx2 = sinTheta2 * std::cos(phi2);
211  G4double diry2 = sinTheta2 * std::sin(phi2);
212  G4double dirz2 = cosTheta2;
213 
214  G4ThreeVector photon1Direction (dirx1,diry1,dirz1);
215  photon1Direction.rotateUz(positronDirection);   
216  // create G4DynamicParticle object for the particle1 
217  G4DynamicParticle* aParticle1= new G4DynamicParticle (G4Gamma::Gamma(),
218                                                           photon1Direction, 
219                                                           photon1Energy);
220  fvect->push_back(aParticle1);
221 
222  G4ThreeVector photon2Direction(dirx2,diry2,dirz2);
223  photon2Direction.rotateUz(positronDirection); 
224     // create G4DynamicParticle object for the particle2
225  G4DynamicParticle* aParticle2= new G4DynamicParticle (G4Gamma::Gamma(),
226                                                           photon2Direction,
227                                                           photon2Energy);
228  fvect->push_back(aParticle2);
229  fParticleChange->SetProposedKineticEnergy(0.);
230  fParticleChange->ProposeTrackStatus(fStopAndKill);
231
232  if (verboseLevel > 1)
233    {
234      G4cout << "-----------------------------------------------------------" << G4endl;
235      G4cout << "Energy balance from G4PenelopeAnnihilation" << G4endl;
236      G4cout << "Kinetic positron energy: " << kineticEnergy/keV << " keV" << G4endl;
237      G4cout << "Total available energy: " << totalAvailableEnergy/keV << " keV " << G4endl;
238      G4cout << "-----------------------------------------------------------" << G4endl;
239      G4cout << "Photon energy 1: " << photon1Energy/keV << " keV" << G4endl;
240      G4cout << "Photon energy 2: " << photon2Energy/keV << " keV" << G4endl;
241      G4cout << "Total final state: " << (photon1Energy+photon2Energy)/keV << 
242        " keV" << G4endl;
243      G4cout << "-----------------------------------------------------------" << G4endl;
244    }
245  if (verboseLevel > 0)
246    {
247     
248      G4double energyDiff = std::fabs(totalAvailableEnergy-photon1Energy-photon2Energy);
249      if (energyDiff > 0.05*keV)
250        G4cout << "Warning from G4PenelopeAnnihilation: problem with energy conservation: " << 
251          (photon1Energy+photon2Energy)/keV << 
252          " keV (final) vs. " << 
253          totalAvailableEnergy/keV << " keV (initial)" << G4endl;
254    }
255  return;
256}
257
258//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
259
260G4double G4PenelopeAnnihilationModel:: ComputeCrossSectionPerElectron(G4double energy)
261{
262  //
263  // Penelope model to calculate cross section for positron annihilation.
264  // The annihilation cross section per electron is calculated according
265  // to the Heitler formula
266  //  W. Heitler, The quantum theory of radiation, Oxford University Press (1954)
267  // in the assumptions of electrons free and at rest.
268  //
269  G4double gamma = 1.0+std::max(energy,1.0*eV)/electron_mass_c2;
270  G4double gamma2 = gamma*gamma;
271  G4double f2 = gamma2-1.0;
272  G4double f1 = std::sqrt(f2);
273  G4double crossSection = fPielr2*((gamma2+4.0*gamma+1.0)*std::log(gamma+f1)/f2
274                         - (gamma+3.0)/f1)/(gamma+1.0);
275  return crossSection;
276}
Note: See TracBrowser for help on using the repository browser.