Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HeatedKleinNishinaCompton.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: G4HeatedKleinNishinaCompton.cc 91726 2015-08-03 15:41:36Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4HeatedKleinNishinaCompton
34 //
35 // Author: Vladimir Grichine on base of M. Maire and V. Ivanchenko code
36 //
37 // Creation date: 15.03.2009
38 //
39 // Modifications: 07.07.2014 V.Ivanchenko make direct inheritence from
40 // G4KleinNishinaCompton
41 //
42 //
43 // Class Description:
44 //
45 // -------------------------------------------------------------------
46 //
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49 
51 #include "globals.hh"
52 #include "G4PhysicalConstants.hh"
53 #include "G4SystemOfUnits.hh"
54 #include "G4RandomDirection.hh"
55 #include "Randomize.hh"
56 #include "G4Log.hh"
57 #include "G4Exp.hh"
58 
59 #include "G4Electron.hh"
60 #include "G4Gamma.hh"
62 
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64 
65 using namespace std;
66 
68  const G4ParticleDefinition* p, const G4String& nam)
69  : G4KleinNishinaCompton(p, nam)
70 {
71  fTemperature = 1.0*keV;
72 }
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 
77 {}
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80 
82  std::vector<G4DynamicParticle*>* fvect,
83  const G4MaterialCutsCouple*,
84  const G4DynamicParticle* aDynamicGamma,
86 {
87  // do nothing below the threshold
88  if(aDynamicGamma->GetKineticEnergy() <= LowEnergyLimit()) { return; }
89 
90  // The scattered gamma energy is sampled according to Klein - Nishina formula.
91  // The random number techniques of Butcher & Messel are used
92  // (Nuc Phys 20(1960),15).
93  // Note : Effects due to binding of atomic electrons are negliged.
94 
95  // We start to prepare a heated electron from Maxwell distribution.
96  // Then we try to boost to the electron rest frame and make scattering.
97  // The final step is to recover new gamma 4momentum in the lab frame
98 
99  G4double eMomentumC2 = G4RandGamma::shoot(1.5, 1.);
100  eMomentumC2 *= 2*electron_mass_c2*fTemperature; // electron (pc)^2
101  G4ThreeVector eMomDir = G4RandomDirection();
102  eMomDir *= std::sqrt(eMomentumC2);
103  G4double eEnergy = std::sqrt(eMomentumC2+electron_mass_c2*electron_mass_c2);
104  G4LorentzVector electron4v = G4LorentzVector(eMomDir,eEnergy);
105  G4ThreeVector bst = electron4v.boostVector();
106 
107  G4LorentzVector gamma4v = aDynamicGamma->Get4Momentum();
108  gamma4v.boost(-bst);
109 
110  G4ThreeVector gammaMomV = gamma4v.vect();
111  G4double gamEnergy0 = gammaMomV.mag();
112 
113 
114  // G4double gamEnergy0 = aDynamicGamma->GetKineticEnergy();
115  G4double E0_m = gamEnergy0 / electron_mass_c2 ;
116 
117  // G4ThreeVector gamDirection0 = /aDynamicGamma->GetMomentumDirection();
118 
119  G4ThreeVector gamDirection0 = gammaMomV/gamEnergy0;
120 
121  // sample the energy rate of the scattered gamma in the electron rest frame
122  //
123 
124  G4double epsilon, epsilonsq, onecost, sint2, greject ;
125 
126  G4double eps0 = 1./(1. + 2.*E0_m);
127  G4double epsilon0sq = eps0*eps0;
128  G4double alpha1 = - G4Log(eps0);
129  G4double alpha2 = 0.5*(1.- epsilon0sq);
130 
131  G4int nloop = 0;
132  do
133  {
134  ++nloop;
135  // false interaction if too many iterations
136  if(nloop > 1000) { return; }
137 
138  if ( alpha1/(alpha1+alpha2) > G4UniformRand() )
139  {
140  epsilon = G4Exp(-alpha1*G4UniformRand()); // eps0**r
141  epsilonsq = epsilon*epsilon;
142 
143  }
144  else
145  {
146  epsilonsq = epsilon0sq + (1.- epsilon0sq)*G4UniformRand();
147  epsilon = sqrt(epsilonsq);
148  };
149 
150  onecost = (1.- epsilon)/(epsilon*E0_m);
151  sint2 = onecost*(2.-onecost);
152  greject = 1. - epsilon*sint2/(1.+ epsilonsq);
153 
154  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
155  } while (greject < G4UniformRand());
156 
157  //
158  // scattered gamma angles. ( Z - axis along the parent gamma)
159  //
160 
161  G4double cosTeta = 1. - onecost;
162  G4double sinTeta = sqrt (sint2);
163  G4double Phi = twopi * G4UniformRand();
164  G4double dirx = sinTeta*cos(Phi), diry = sinTeta*sin(Phi), dirz = cosTeta;
165 
166  //
167  // update G4VParticleChange for the scattered gamma
168  //
169 
170  G4ThreeVector gamDirection1 ( dirx,diry,dirz );
171  gamDirection1.rotateUz(gamDirection0);
172  G4double gamEnergy1 = epsilon*gamEnergy0;
173  gamDirection1 *= gamEnergy1;
174 
175  G4LorentzVector gamma4vfinal = G4LorentzVector(gamDirection1,gamEnergy1);
176 
177 
178  // kinematic of the scattered electron
179  //
180 
181  G4double eKinEnergy = gamEnergy0 - gamEnergy1;
182  G4ThreeVector eDirection = gamEnergy0*gamDirection0 - gamEnergy1*gamDirection1;
183  eDirection = eDirection.unit();
184  G4double eFinalMom = std::sqrt(eKinEnergy*(eKinEnergy+2*electron_mass_c2));
185  eDirection *= eFinalMom;
186  G4LorentzVector e4vfinal = G4LorentzVector(eDirection,gamEnergy1+electron_mass_c2);
187 
188  gamma4vfinal.boost(bst);
189  e4vfinal.boost(bst);
190 
191  gamDirection1 = gamma4vfinal.vect();
192  gamEnergy1 = gamDirection1.mag();
193  gamDirection1 /= gamEnergy1;
194 
195  G4double edep = 0.0;
196  if(gamEnergy1 > lowestSecondaryEnergy) {
199  } else {
202  edep = gamEnergy1;
203  }
204 
205  //
206  // kinematic of the scattered electron
207  //
208  eKinEnergy = e4vfinal.t()-electron_mass_c2;
209 
210  if(eKinEnergy > lowestSecondaryEnergy) {
211 
212  eDirection = e4vfinal.vect().unit();
213 
214  // create G4DynamicParticle object for the electron.
215  G4DynamicParticle* dp =
216  new G4DynamicParticle(theElectron,eDirection,eKinEnergy);
217  fvect->push_back(dp);
218  } else {
219  edep += eKinEnergy;
220  }
221  // energy balance
222  if(edep > 0.0) {
224  }
225 }
226 
227 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
228 
229 
ThreeVector shoot(const G4int Ap, const G4int Af)
Hep3Vector boostVector() const
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
G4ParticleChangeForGamma * fParticleChange
G4double GetKineticEnergy() const
const char * p
Definition: xmltok.h:285
G4ThreeVector G4RandomDirection()
int G4int
Definition: G4Types.hh:78
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double twopi
Definition: G4SIunits.hh:76
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
HepLorentzVector & boost(double, double, double)
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) final
float electron_mass_c2
Definition: hepunit.py:274
G4LorentzVector Get4Momentum() const
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
Hep3Vector unit() const
G4HeatedKleinNishinaCompton(const G4ParticleDefinition *p=nullptr, const G4String &nam="Heated-Klein-Nishina")
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4ParticleDefinition * theElectron
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
double mag() const
static constexpr double keV
Definition: G4SIunits.hh:216
double epsilon(double density, double temperature)
CLHEP::HepLorentzVector G4LorentzVector