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