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 ]

Diff markup

Differences between /processes/electromagnetic/lowenergy/src/G4PenelopeAnnihilationModel.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4PenelopeAnnihilationModel.cc (Version 9.5.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 // $Id: G4PenelopeAnnihilationModel.cc,v 1.4 2009-06-10 13:32:36 mantero Exp $
                                                   >>  27 // GEANT4 tag $Name: not supported by cvs2svn $
 26 //                                                 28 //
 27 // Author: Luciano Pandola                         29 // Author: Luciano Pandola
 28 //                                                 30 //
 29 // History:                                        31 // History:
 30 // --------                                        32 // --------
 31 // 29 Oct 2008   L Pandola    Migration from p     33 // 29 Oct 2008   L Pandola    Migration from process to model 
 32 // 15 Apr 2009   V Ivanchenko Cleanup initiali <<  34 // 15 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
 33 //                    secondaries:             << 
 34 //                  - apply internal high-ener     35 //                  - apply internal high-energy limit only in constructor 
 35 //                  - do not apply low-energy      36 //                  - do not apply low-energy limit (default is 0)
 36 //                  - do not use G4ElementSele     37 //                  - do not use G4ElementSelector
 37 // 02 Oct 2013   L.Pandola    Migration to MT  << 
 38                                                    38 
 39 #include "G4PenelopeAnnihilationModel.hh"          39 #include "G4PenelopeAnnihilationModel.hh"
 40 #include "G4PhysicalConstants.hh"              << 
 41 #include "G4SystemOfUnits.hh"                  << 
 42 #include "G4ParticleDefinition.hh"                 40 #include "G4ParticleDefinition.hh"
 43 #include "G4MaterialCutsCouple.hh"                 41 #include "G4MaterialCutsCouple.hh"
 44 #include "G4ProductionCutsTable.hh"                42 #include "G4ProductionCutsTable.hh"
 45 #include "G4DynamicParticle.hh"                    43 #include "G4DynamicParticle.hh"
 46 #include "G4Gamma.hh"                              44 #include "G4Gamma.hh"
 47                                                    45 
 48 //....oooOO0OOooo........oooOO0OOooo........oo     46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 49                                                    47 
 50 G4double G4PenelopeAnnihilationModel::fPielr2  << 
 51                                                    48 
 52 G4PenelopeAnnihilationModel::G4PenelopeAnnihil <<  49 G4PenelopeAnnihilationModel::G4PenelopeAnnihilationModel(const G4ParticleDefinition*,
 53                                              c     50                                              const G4String& nam)
 54   :G4VEmModel(nam),fParticleChange(nullptr),fP <<  51   :G4VEmModel(nam),fParticleChange(0),isInitialised(false)
 55 {                                                  52 {
 56   fIntrinsicLowEnergyLimit = 0.0;                  53   fIntrinsicLowEnergyLimit = 0.0;
 57   fIntrinsicHighEnergyLimit = 100.0*GeV;           54   fIntrinsicHighEnergyLimit = 100.0*GeV;
                                                   >>  55   //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
 58   SetHighEnergyLimit(fIntrinsicHighEnergyLimit     56   SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
 59                                                <<  57  
 60   if (part)                                    << 
 61     SetParticle(part);                         << 
 62                                                << 
 63   //Calculate variable that will be used later     58   //Calculate variable that will be used later on
 64   fPielr2 = pi*classic_electr_radius*classic_e     59   fPielr2 = pi*classic_electr_radius*classic_electr_radius;
 65                                                    60 
 66   fVerboseLevel= 0;                            <<  61   verboseLevel= 0;
 67   // Verbosity scale:                              62   // Verbosity scale:
 68   // 0 = nothing                                   63   // 0 = nothing 
 69   // 1 = warning for energy non-conservation       64   // 1 = warning for energy non-conservation 
 70   // 2 = details of energy budget                  65   // 2 = details of energy budget
 71   // 3 = calculation of cross sections, file o     66   // 3 = calculation of cross sections, file openings, sampling of atoms
 72   // 4 = entering in methods                       67   // 4 = entering in methods
                                                   >>  68 
 73 }                                                  69 }
 74                                                    70 
 75 //....oooOO0OOooo........oooOO0OOooo........oo     71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 76                                                    72 
 77 G4PenelopeAnnihilationModel::~G4PenelopeAnnihi     73 G4PenelopeAnnihilationModel::~G4PenelopeAnnihilationModel()
 78 {;}                                                74 {;}
 79                                                    75 
 80 //....oooOO0OOooo........oooOO0OOooo........oo     76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 81                                                    77 
 82 void G4PenelopeAnnihilationModel::Initialise(c <<  78 void G4PenelopeAnnihilationModel::Initialise(const G4ParticleDefinition*,
 83                const G4DataVector&)                79                const G4DataVector&)
 84 {                                                  80 {
 85   if (fVerboseLevel > 3)                       <<  81   if (verboseLevel > 3)
 86     G4cout << "Calling G4PenelopeAnnihilationM     82     G4cout << "Calling G4PenelopeAnnihilationModel::Initialise()" << G4endl;
 87   SetParticle(part);                           << 
 88                                                << 
 89   if (IsMaster() && part == fParticle)         << 
 90     {                                          << 
 91                                                    83 
 92       if(fVerboseLevel > 0) {                  <<  84   if(verboseLevel > 0) {
 93   G4cout << "Penelope Annihilation model is in <<  85     G4cout << "Penelope Annihilation model is initialized " << G4endl
 94          << "Energy range: "                   <<  86      << "Energy range: "
 95          << LowEnergyLimit() / keV << " keV -  <<  87      << LowEnergyLimit() / keV << " keV - "
 96          << HighEnergyLimit() / GeV << " GeV"  <<  88      << HighEnergyLimit() / GeV << " GeV"
 97          << G4endl;                            <<  89      << G4endl;
 98       }                                        <<  90   }
 99     }                                          << 
100                                                    91 
101   if(fIsInitialised) return;                   <<  92   if(isInitialised) return;
102   fParticleChange = GetParticleChangeForGamma(     93   fParticleChange = GetParticleChangeForGamma();
103   fIsInitialised = true;                       <<  94   isInitialised = true; 
104 }                                              << 
105                                                << 
106 //....oooOO0OOooo........oooOO0OOooo........oo << 
107 void G4PenelopeAnnihilationModel::InitialiseLo << 
108               G4VEmModel* masterModel)         << 
109 {                                              << 
110   if (fVerboseLevel > 3)                       << 
111     G4cout << "Calling G4PenelopeAnnihilationM << 
112                                                << 
113   //                                           << 
114   //Check that particle matches: one might hav << 
115   //for e+ and e-).                            << 
116   //                                           << 
117   if (part == fParticle)                       << 
118     {                                          << 
119       //Get the const table pointers from the  << 
120       const G4PenelopeAnnihilationModel* theMo << 
121         static_cast<G4PenelopeAnnihilationMode << 
122                                                << 
123       //Same verbosity for all workers, as the << 
124       fVerboseLevel = theModel->fVerboseLevel; << 
125     }                                          << 
126 }                                                  95 }
127                                                    96 
128 //....oooOO0OOooo........oooOO0OOooo........oo     97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
129                                                    98 
130 G4double G4PenelopeAnnihilationModel::ComputeC     99 G4double G4PenelopeAnnihilationModel::ComputeCrossSectionPerAtom(
131                                        const G    100                                        const G4ParticleDefinition*,
132                                              G    101                                              G4double energy,
133                                              G    102                                              G4double Z, G4double,
134                                              G    103                                              G4double, G4double)
135 {                                                 104 {
136   if (fVerboseLevel > 3)                       << 105   if (verboseLevel > 3)
137     G4cout << "Calling ComputeCrossSectionPerA    106     G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopeAnnihilationModel" << 
138       G4endl;                                     107       G4endl;
139                                                   108 
140   G4double cs = Z*ComputeCrossSectionPerElectr    109   G4double cs = Z*ComputeCrossSectionPerElectron(energy);
141                                                   110   
142   if (fVerboseLevel > 2)                       << 111   if (verboseLevel > 2)
143     G4cout << "Annihilation cross Section at "    112     G4cout << "Annihilation cross Section at " << energy/keV << " keV for Z=" << Z << 
144       " = " << cs/barn << " barn" << G4endl;      113       " = " << cs/barn << " barn" << G4endl;
145   return cs;                                      114   return cs;
146 }                                                 115 }
147                                                   116 
148 //....oooOO0OOooo........oooOO0OOooo........oo    117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
149                                                   118 
150 void G4PenelopeAnnihilationModel::SampleSecond    119 void G4PenelopeAnnihilationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
151                 const G4MaterialCutsCouple*,      120                 const G4MaterialCutsCouple*,
152                 const G4DynamicParticle* aDyna    121                 const G4DynamicParticle* aDynamicPositron,
153                 G4double,                         122                 G4double,
154                 G4double)                         123                 G4double)
155 {                                                 124 {
156   //                                              125   //
157   // Penelope model to sample final state for     126   // Penelope model to sample final state for positron annihilation. 
158   // Target eletrons are assumed to be free an    127   // Target eletrons are assumed to be free and at rest. Binding effects enabling 
159   // one-photon annihilation are neglected.       128   // one-photon annihilation are neglected.
160   // For annihilation at rest, two back-to-bac    129   // For annihilation at rest, two back-to-back photons are emitted, having energy of 511 keV 
161   // and isotropic angular distribution.          130   // and isotropic angular distribution.
162   // For annihilation in flight, it is used th    131   // For annihilation in flight, it is used the theory from 
163   //  W. Heitler, The quantum theory of radiat    132   //  W. Heitler, The quantum theory of radiation, Oxford University Press (1954)
164   // The two photons can have different energy    133   // The two photons can have different energy. The efficiency of the sampling algorithm 
165   // of the photon energy from the dSigma/dE d    134   // of the photon energy from the dSigma/dE distribution is practically 100% for 
166   // positrons of kinetic energy < 10 keV. It     135   // positrons of kinetic energy < 10 keV. It reaches a minimum (about 80%) at energy 
167   // of about 10 MeV.                             136   // of about 10 MeV.
168   // The angle theta is kinematically linked t    137   // The angle theta is kinematically linked to the photon energy, to ensure momentum 
169   // conservation. The angle phi is sampled is    138   // conservation. The angle phi is sampled isotropically for the first gamma.
170   //                                              139   //
171   if (fVerboseLevel > 3)                       << 140   if (verboseLevel > 3)
172     G4cout << "Calling SamplingSecondaries() o    141     G4cout << "Calling SamplingSecondaries() of G4PenelopeAnnihilationModel" << G4endl;
173                                                   142 
174   G4double kineticEnergy = aDynamicPositron->G    143   G4double kineticEnergy = aDynamicPositron->GetKineticEnergy();
175                                                   144 
176   // kill primary                                 145   // kill primary
177   fParticleChange->SetProposedKineticEnergy(0.    146   fParticleChange->SetProposedKineticEnergy(0.);
178   fParticleChange->ProposeTrackStatus(fStopAnd    147   fParticleChange->ProposeTrackStatus(fStopAndKill);
179                                                   148   
180   if (kineticEnergy == 0.0)                       149   if (kineticEnergy == 0.0)
181     {                                             150     {
182       //Old AtRestDoIt                            151       //Old AtRestDoIt
183       G4double cosTheta = -1.0+2.0*G4UniformRa    152       G4double cosTheta = -1.0+2.0*G4UniformRand();
184       G4double sinTheta = std::sqrt(1.0-cosThe    153       G4double sinTheta = std::sqrt(1.0-cosTheta*cosTheta);
185       G4double phi = twopi*G4UniformRand();       154       G4double phi = twopi*G4UniformRand();
186       G4ThreeVector direction (sinTheta*std::c    155       G4ThreeVector direction (sinTheta*std::cos(phi),sinTheta*std::sin(phi),cosTheta);
187       G4DynamicParticle* firstGamma = new G4Dy    156       G4DynamicParticle* firstGamma = new G4DynamicParticle (G4Gamma::Gamma(),
188                    direction, electron_mass_c2    157                    direction, electron_mass_c2);
189       G4DynamicParticle* secondGamma = new G4D    158       G4DynamicParticle* secondGamma = new G4DynamicParticle (G4Gamma::Gamma(),
190                     -direction, electron_mass_    159                     -direction, electron_mass_c2);
191                                                   160   
192       fvect->push_back(firstGamma);               161       fvect->push_back(firstGamma);
193       fvect->push_back(secondGamma);              162       fvect->push_back(secondGamma);
194       return;                                     163       return;
195     }                                             164     }
196                                                   165 
197   //This is the "PostStep" case (annihilation     166   //This is the "PostStep" case (annihilation in flight)
198   G4ParticleMomentum positronDirection =          167   G4ParticleMomentum positronDirection = 
199     aDynamicPositron->GetMomentumDirection();     168     aDynamicPositron->GetMomentumDirection();
200   G4double gamma = 1.0 + std::max(kineticEnerg    169   G4double gamma = 1.0 + std::max(kineticEnergy,1.0*eV)/electron_mass_c2;
201   G4double gamma21 = std::sqrt(gamma*gamma-1);    170   G4double gamma21 = std::sqrt(gamma*gamma-1);
202   G4double ani = 1.0+gamma;                       171   G4double ani = 1.0+gamma;
203   G4double chimin = 1.0/(ani+gamma21);            172   G4double chimin = 1.0/(ani+gamma21);
204   G4double rchi = (1.0-chimin)/chimin;            173   G4double rchi = (1.0-chimin)/chimin;
205   G4double gt0 = ani*ani-2.0;                     174   G4double gt0 = ani*ani-2.0;
206   G4double test=0.0;                              175   G4double test=0.0;
207   G4double epsilon = 0;                           176   G4double epsilon = 0;
208   do{                                             177   do{
209     epsilon = chimin*std::pow(rchi,G4UniformRa    178     epsilon = chimin*std::pow(rchi,G4UniformRand());
210     G4double reject = ani*ani*(1.0-epsilon)+2.    179     G4double reject = ani*ani*(1.0-epsilon)+2.0*gamma-(1.0/epsilon);
211     test = G4UniformRand()*gt0-reject;            180     test = G4UniformRand()*gt0-reject;
212   }while(test>0);                                 181   }while(test>0);
213                                                   182    
214   G4double totalAvailableEnergy = kineticEnerg    183   G4double totalAvailableEnergy = kineticEnergy + 2.0*electron_mass_c2;
215   G4double photon1Energy = epsilon*totalAvaila    184   G4double photon1Energy = epsilon*totalAvailableEnergy;
216   G4double photon2Energy = (1.0-epsilon)*total    185   G4double photon2Energy = (1.0-epsilon)*totalAvailableEnergy;
217   G4double cosTheta1 = (ani-1.0/epsilon)/gamma    186   G4double cosTheta1 = (ani-1.0/epsilon)/gamma21;
218   G4double cosTheta2 = (ani-1.0/(1.0-epsilon))    187   G4double cosTheta2 = (ani-1.0/(1.0-epsilon))/gamma21;
219                                                   188   
                                                   >> 189   //G4double localEnergyDeposit = 0.; 
                                                   >> 190 
220   G4double sinTheta1 = std::sqrt(1.-cosTheta1*    191   G4double sinTheta1 = std::sqrt(1.-cosTheta1*cosTheta1);
221   G4double phi1  = twopi * G4UniformRand();       192   G4double phi1  = twopi * G4UniformRand();
222   G4double dirx1 = sinTheta1 * std::cos(phi1);    193   G4double dirx1 = sinTheta1 * std::cos(phi1);
223   G4double diry1 = sinTheta1 * std::sin(phi1);    194   G4double diry1 = sinTheta1 * std::sin(phi1);
224   G4double dirz1 = cosTheta1;                     195   G4double dirz1 = cosTheta1;
225                                                   196   
226   G4double sinTheta2 = std::sqrt(1.-cosTheta2*    197   G4double sinTheta2 = std::sqrt(1.-cosTheta2*cosTheta2);
227   G4double phi2  = phi1+pi;                       198   G4double phi2  = phi1+pi;
228   G4double dirx2 = sinTheta2 * std::cos(phi2);    199   G4double dirx2 = sinTheta2 * std::cos(phi2);
229   G4double diry2 = sinTheta2 * std::sin(phi2);    200   G4double diry2 = sinTheta2 * std::sin(phi2);
230   G4double dirz2 = cosTheta2;                     201   G4double dirz2 = cosTheta2;
231                                                   202   
232   G4ThreeVector photon1Direction (dirx1,diry1,    203   G4ThreeVector photon1Direction (dirx1,diry1,dirz1);
233   photon1Direction.rotateUz(positronDirection)    204   photon1Direction.rotateUz(positronDirection);   
234   // create G4DynamicParticle object for the p    205   // create G4DynamicParticle object for the particle1  
235   G4DynamicParticle* aParticle1= new G4Dynamic    206   G4DynamicParticle* aParticle1= new G4DynamicParticle (G4Gamma::Gamma(),
236                  photon1Direction,                207                  photon1Direction, 
237                  photon1Energy);                  208                  photon1Energy);
238   fvect->push_back(aParticle1);                   209   fvect->push_back(aParticle1);
239                                                   210  
240   G4ThreeVector photon2Direction(dirx2,diry2,d    211   G4ThreeVector photon2Direction(dirx2,diry2,dirz2);
241   photon2Direction.rotateUz(positronDirection)    212   photon2Direction.rotateUz(positronDirection); 
242   // create G4DynamicParticle object for the p << 213      // create G4DynamicParticle object for the particle2 
243   G4DynamicParticle* aParticle2= new G4Dynamic    214   G4DynamicParticle* aParticle2= new G4DynamicParticle (G4Gamma::Gamma(),
244                  photon2Direction,                215                  photon2Direction,
245                  photon2Energy);                  216                  photon2Energy);
246   fvect->push_back(aParticle2);                   217   fvect->push_back(aParticle2);
247                                                   218 
248   if (fVerboseLevel > 1)                       << 219   if (verboseLevel > 1)
249     {                                             220     {
250       G4cout << "-----------------------------    221       G4cout << "-----------------------------------------------------------" << G4endl;
251       G4cout << "Energy balance from G4Penelop    222       G4cout << "Energy balance from G4PenelopeAnnihilation" << G4endl;
252       G4cout << "Kinetic positron energy: " <<    223       G4cout << "Kinetic positron energy: " << kineticEnergy/keV << " keV" << G4endl;
253       G4cout << "Total available energy: " <<     224       G4cout << "Total available energy: " << totalAvailableEnergy/keV << " keV " << G4endl;
254       G4cout << "-----------------------------    225       G4cout << "-----------------------------------------------------------" << G4endl;
255       G4cout << "Photon energy 1: " << photon1    226       G4cout << "Photon energy 1: " << photon1Energy/keV << " keV" << G4endl;
256       G4cout << "Photon energy 2: " << photon2    227       G4cout << "Photon energy 2: " << photon2Energy/keV << " keV" << G4endl;
257       G4cout << "Total final state: " << (phot    228       G4cout << "Total final state: " << (photon1Energy+photon2Energy)/keV << 
258   " keV" << G4endl;                               229   " keV" << G4endl;
259       G4cout << "-----------------------------    230       G4cout << "-----------------------------------------------------------" << G4endl;
260     }                                             231     }
261   if (fVerboseLevel > 0)                       << 232   if (verboseLevel > 0)
262     {                                             233     {      
263       G4double energyDiff = std::fabs(totalAva    234       G4double energyDiff = std::fabs(totalAvailableEnergy-photon1Energy-photon2Energy);
264       if (energyDiff > 0.05*keV)                  235       if (energyDiff > 0.05*keV)
265   G4cout << "Warning from G4PenelopeAnnihilati    236   G4cout << "Warning from G4PenelopeAnnihilation: problem with energy conservation: " << 
266     (photon1Energy+photon2Energy)/keV <<          237     (photon1Energy+photon2Energy)/keV << 
267     " keV (final) vs. " <<                        238     " keV (final) vs. " << 
268     totalAvailableEnergy/keV << " keV (initial    239     totalAvailableEnergy/keV << " keV (initial)" << G4endl;
269     }                                             240     }
270   return;                                         241   return;
271 }                                                 242 }
272                                                   243 
273 //....oooOO0OOooo........oooOO0OOooo........oo    244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
274                                                   245 
275 G4double G4PenelopeAnnihilationModel:: Compute    246 G4double G4PenelopeAnnihilationModel:: ComputeCrossSectionPerElectron(G4double energy)
276 {                                                 247 {
277   //                                              248   //
278   // Penelope model to calculate cross section    249   // Penelope model to calculate cross section for positron annihilation.
279   // The annihilation cross section per electr    250   // The annihilation cross section per electron is calculated according 
280   // to the Heitler formula                       251   // to the Heitler formula
281   //  W. Heitler, The quantum theory of radiat    252   //  W. Heitler, The quantum theory of radiation, Oxford University Press (1954)
282   // in the assumptions of electrons free and     253   // in the assumptions of electrons free and at rest.
283   //                                              254   //
284   G4double gamma = 1.0+std::max(energy,1.0*eV)    255   G4double gamma = 1.0+std::max(energy,1.0*eV)/electron_mass_c2;
285   G4double gamma2 = gamma*gamma;                  256   G4double gamma2 = gamma*gamma;
286   G4double f2 = gamma2-1.0;                       257   G4double f2 = gamma2-1.0;
287   G4double f1 = std::sqrt(f2);                    258   G4double f1 = std::sqrt(f2);
288   G4double crossSection = fPielr2*((gamma2+4.0 << 259   G4double crossSection = fPielr2*((gamma2+4.0*gamma+1.0)*std::log(gamma+f1)/f2
289        - (gamma+3.0)/f1)/(gamma+1.0);             260        - (gamma+3.0)/f1)/(gamma+1.0);
290   return crossSection;                            261   return crossSection;
291 }                                              << 
292                                                << 
293 //....oooOO0OOooo........oooOO0OOooo........oo << 
294                                                << 
295 void G4PenelopeAnnihilationModel::SetParticle( << 
296 {                                              << 
297   if(!fParticle) {                             << 
298     fParticle = p;                             << 
299   }                                            << 
300 }                                                 262 }
301                                                   263