Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4PenelopeAnnihilationModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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