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