Geant4  10.01.p01
G4PenelopeAnnihilationModel.cc
Go to the documentation of this file.
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 74452 2013-10-07 15:08:00Z gcosmo $
27 //
28 // Author: Luciano Pandola
29 //
30 // History:
31 // --------
32 // 29 Oct 2008 L Pandola Migration from process to model
33 // 15 Apr 2009 V Ivanchenko Cleanup initialisation and generation of
34 // secondaries:
35 // - apply internal high-energy limit only in constructor
36 // - do not apply low-energy limit (default is 0)
37 // - do not use G4ElementSelector
38 // 02 Oct 2013 L.Pandola Migration to MT
39 
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4ParticleDefinition.hh"
44 #include "G4MaterialCutsCouple.hh"
45 #include "G4ProductionCutsTable.hh"
46 #include "G4DynamicParticle.hh"
47 #include "G4Gamma.hh"
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50 
52 
54  const G4String& nam)
55  :G4VEmModel(nam),fParticleChange(0),fParticle(0),isInitialised(false)
56 {
59  // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
61 
62  if (part)
63  SetParticle(part);
64 
65  //Calculate variable that will be used later on
66  fPielr2 = pi*classic_electr_radius*classic_electr_radius;
67 
68  verboseLevel= 0;
69  // Verbosity scale:
70  // 0 = nothing
71  // 1 = warning for energy non-conservation
72  // 2 = details of energy budget
73  // 3 = calculation of cross sections, file openings, sampling of atoms
74  // 4 = entering in methods
75 
76 }
77 
78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
79 
81 {;}
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84 
86  const G4DataVector&)
87 {
88  if (verboseLevel > 3)
89  G4cout << "Calling G4PenelopeAnnihilationModel::Initialise()" << G4endl;
90  SetParticle(part);
91 
92  if (IsMaster() && part == fParticle)
93  {
94 
95  if(verboseLevel > 0) {
96  G4cout << "Penelope Annihilation model is initialized " << G4endl
97  << "Energy range: "
98  << LowEnergyLimit() / keV << " keV - "
99  << HighEnergyLimit() / GeV << " GeV"
100  << G4endl;
101  }
102  }
103 
104  if(isInitialised) return;
106  isInitialised = true;
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111  G4VEmModel* masterModel)
112 {
113  if (verboseLevel > 3)
114  G4cout << "Calling G4PenelopeAnnihilationModel::InitialiseLocal()" << G4endl;
115 
116  //
117  //Check that particle matches: one might have multiple master models (e.g.
118  //for e+ and e-).
119  //
120  if (part == fParticle)
121  {
122  //Get the const table pointers from the master to the workers
123  const G4PenelopeAnnihilationModel* theModel =
124  static_cast<G4PenelopeAnnihilationModel*> (masterModel);
125 
126  //Same verbosity for all workers, as the master
127  verboseLevel = theModel->verboseLevel;
128  }
129 }
130 
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
132 
134  const G4ParticleDefinition*,
136  G4double Z, G4double,
138 {
139  if (verboseLevel > 3)
140  G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopeAnnihilationModel" <<
141  G4endl;
142 
144 
145  if (verboseLevel > 2)
146  G4cout << "Annihilation cross Section at " << energy/keV << " keV for Z=" << Z <<
147  " = " << cs/barn << " barn" << G4endl;
148  return cs;
149 }
150 
151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
152 
153 void G4PenelopeAnnihilationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
154  const G4MaterialCutsCouple*,
155  const G4DynamicParticle* aDynamicPositron,
156  G4double,
157  G4double)
158 {
159  //
160  // Penelope model to sample final state for positron annihilation.
161  // Target eletrons are assumed to be free and at rest. Binding effects enabling
162  // one-photon annihilation are neglected.
163  // For annihilation at rest, two back-to-back photons are emitted, having energy of 511 keV
164  // and isotropic angular distribution.
165  // For annihilation in flight, it is used the theory from
166  // W. Heitler, The quantum theory of radiation, Oxford University Press (1954)
167  // The two photons can have different energy. The efficiency of the sampling algorithm
168  // of the photon energy from the dSigma/dE distribution is practically 100% for
169  // positrons of kinetic energy < 10 keV. It reaches a minimum (about 80%) at energy
170  // of about 10 MeV.
171  // The angle theta is kinematically linked to the photon energy, to ensure momentum
172  // conservation. The angle phi is sampled isotropically for the first gamma.
173  //
174  if (verboseLevel > 3)
175  G4cout << "Calling SamplingSecondaries() of G4PenelopeAnnihilationModel" << G4endl;
176 
177  G4double kineticEnergy = aDynamicPositron->GetKineticEnergy();
178 
179  // kill primary
182 
183  if (kineticEnergy == 0.0)
184  {
185  //Old AtRestDoIt
186  G4double cosTheta = -1.0+2.0*G4UniformRand();
187  G4double sinTheta = std::sqrt(1.0-cosTheta*cosTheta);
188  G4double phi = twopi*G4UniformRand();
189  G4ThreeVector direction (sinTheta*std::cos(phi),sinTheta*std::sin(phi),cosTheta);
191  direction, electron_mass_c2);
192  G4DynamicParticle* secondGamma = new G4DynamicParticle (G4Gamma::Gamma(),
193  -direction, electron_mass_c2);
194 
195  fvect->push_back(firstGamma);
196  fvect->push_back(secondGamma);
197  return;
198  }
199 
200  //This is the "PostStep" case (annihilation in flight)
201  G4ParticleMomentum positronDirection =
202  aDynamicPositron->GetMomentumDirection();
203  G4double gamma = 1.0 + std::max(kineticEnergy,1.0*eV)/electron_mass_c2;
204  G4double gamma21 = std::sqrt(gamma*gamma-1);
205  G4double ani = 1.0+gamma;
206  G4double chimin = 1.0/(ani+gamma21);
207  G4double rchi = (1.0-chimin)/chimin;
208  G4double gt0 = ani*ani-2.0;
209  G4double test=0.0;
210  G4double epsilon = 0;
211  do{
212  epsilon = chimin*std::pow(rchi,G4UniformRand());
213  G4double reject = ani*ani*(1.0-epsilon)+2.0*gamma-(1.0/epsilon);
214  test = G4UniformRand()*gt0-reject;
215  }while(test>0);
216 
217  G4double totalAvailableEnergy = kineticEnergy + 2.0*electron_mass_c2;
218  G4double photon1Energy = epsilon*totalAvailableEnergy;
219  G4double photon2Energy = (1.0-epsilon)*totalAvailableEnergy;
220  G4double cosTheta1 = (ani-1.0/epsilon)/gamma21;
221  G4double cosTheta2 = (ani-1.0/(1.0-epsilon))/gamma21;
222 
223  //G4double localEnergyDeposit = 0.;
224 
225  G4double sinTheta1 = std::sqrt(1.-cosTheta1*cosTheta1);
226  G4double phi1 = twopi * G4UniformRand();
227  G4double dirx1 = sinTheta1 * std::cos(phi1);
228  G4double diry1 = sinTheta1 * std::sin(phi1);
229  G4double dirz1 = cosTheta1;
230 
231  G4double sinTheta2 = std::sqrt(1.-cosTheta2*cosTheta2);
232  G4double phi2 = phi1+pi;
233  G4double dirx2 = sinTheta2 * std::cos(phi2);
234  G4double diry2 = sinTheta2 * std::sin(phi2);
235  G4double dirz2 = cosTheta2;
236 
237  G4ThreeVector photon1Direction (dirx1,diry1,dirz1);
238  photon1Direction.rotateUz(positronDirection);
239  // create G4DynamicParticle object for the particle1
241  photon1Direction,
242  photon1Energy);
243  fvect->push_back(aParticle1);
244 
245  G4ThreeVector photon2Direction(dirx2,diry2,dirz2);
246  photon2Direction.rotateUz(positronDirection);
247  // create G4DynamicParticle object for the particle2
249  photon2Direction,
250  photon2Energy);
251  fvect->push_back(aParticle2);
252 
253  if (verboseLevel > 1)
254  {
255  G4cout << "-----------------------------------------------------------" << G4endl;
256  G4cout << "Energy balance from G4PenelopeAnnihilation" << G4endl;
257  G4cout << "Kinetic positron energy: " << kineticEnergy/keV << " keV" << G4endl;
258  G4cout << "Total available energy: " << totalAvailableEnergy/keV << " keV " << G4endl;
259  G4cout << "-----------------------------------------------------------" << G4endl;
260  G4cout << "Photon energy 1: " << photon1Energy/keV << " keV" << G4endl;
261  G4cout << "Photon energy 2: " << photon2Energy/keV << " keV" << G4endl;
262  G4cout << "Total final state: " << (photon1Energy+photon2Energy)/keV <<
263  " keV" << G4endl;
264  G4cout << "-----------------------------------------------------------" << G4endl;
265  }
266  if (verboseLevel > 0)
267  {
268  G4double energyDiff = std::fabs(totalAvailableEnergy-photon1Energy-photon2Energy);
269  if (energyDiff > 0.05*keV)
270  G4cout << "Warning from G4PenelopeAnnihilation: problem with energy conservation: " <<
271  (photon1Energy+photon2Energy)/keV <<
272  " keV (final) vs. " <<
273  totalAvailableEnergy/keV << " keV (initial)" << G4endl;
274  }
275  return;
276 }
277 
278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
279 
281 {
282  //
283  // Penelope model to calculate cross section for positron annihilation.
284  // The annihilation cross section per electron is calculated according
285  // to the Heitler formula
286  // W. Heitler, The quantum theory of radiation, Oxford University Press (1954)
287  // in the assumptions of electrons free and at rest.
288  //
289  G4double gamma = 1.0+std::max(energy,1.0*eV)/electron_mass_c2;
290  G4double gamma2 = gamma*gamma;
291  G4double f2 = gamma2-1.0;
292  G4double f1 = std::sqrt(f2);
293  G4double crossSection = fPielr2*((gamma2+4.0*gamma+1.0)*std::log(gamma+f1)/f2
294  - (gamma+3.0)/f1)/(gamma+1.0);
295  return crossSection;
296 }
297 
298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
299 
301 {
302  if(!fParticle) {
303  fParticle = p;
304  }
305 }
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:622
static const G4double f2
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:615
const G4double pi
G4double ComputeCrossSectionPerElectron(G4double energy)
G4bool IsMaster() const
Definition: G4VEmModel.hh:699
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:706
const G4ParticleDefinition * fParticle
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
void SetParticle(const G4ParticleDefinition *)
#define G4UniformRand()
Definition: Randomize.hh:95
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
static const double GeV
Definition: G4SIunits.hh:196
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
static const G4double f1
static const double eV
Definition: G4SIunits.hh:194
G4PenelopeAnnihilationModel(const G4ParticleDefinition *p=0, const G4String &processName="PenAnnih")
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *)
G4ParticleChangeForGamma * fParticleChange
static const double barn
Definition: G4SIunits.hh:95
static const double keV
Definition: G4SIunits.hh:195
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
G4ThreeVector G4ParticleMomentum
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:134