source: trunk/source/processes/electromagnetic/standard/src/G4PolarizedComptonScattering.cc@ 1005

Last change on this file since 1005 was 991, checked in by garnier, 17 years ago

update
File size: 10.4 KB
Line 
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// $Id: G4PolarizedComptonScattering.cc,v 1.18 2008/10/15 17:53:44 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-02 $
29//
30//
31//---------- G4PolarizedComptonScattering physics process ----------------------
32// by Vicente Lara, March 1998
33//
34// -----------------------------------------------------------------------------
35// Corrections by Rui Curado da Silva (Nov. 2000)
36// - Sampling of Phi
37// - Depolarization probability
38//
39// 13-07-01, DoIt: suppression of production cut for the electron (mma)
40// 20-09-01, DoIt: fminimalEnergy = 1*eV (mma)
41// 04-05-05, Inheritance from ComptonScattering52 (V.Ivanchenko)
42// 30-01-06, DoIt : return G4ComptonScattering52::PostStepDoIt(aTrack,aStep) mma
43//
44// -----------------------------------------------------------------------------
45
46#include "G4PolarizedComptonScattering.hh"
47
48//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49
50using namespace std;
51
52G4PolarizedComptonScattering::G4PolarizedComptonScattering(const G4String& pname)
53 : G4ComptonScattering52 (pname)
54{ }
55
56//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57
58G4VParticleChange* G4PolarizedComptonScattering::PostStepDoIt(
59 const G4Track& aTrack,
60 const G4Step& aStep)
61//
62// The scattered gamma energy is sampled according to Klein - Nishina formula.
63// The random number techniques of Butcher & Messel are used
64// (Nuc Phys 20(1960),15).
65// GEANT4 internal units
66//
67// Note : Effects due to binding of atomic electrons are negliged.
68
69{
70 aParticleChange.Initialize(aTrack);
71
72 const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
73
74 G4ThreeVector GammaPolarization0 = aDynamicGamma->GetPolarization();
75
76 if (std::abs(GammaPolarization0.mag() - 1.e0) > 1.e-14)
77 return G4ComptonScattering52::PostStepDoIt(aTrack,aStep);
78
79 G4double GammaEnergy0 = aDynamicGamma->GetKineticEnergy();
80 G4double E0_m = GammaEnergy0 / electron_mass_c2;
81
82 G4ParticleMomentum GammaDirection0 = aDynamicGamma->GetMomentumDirection();
83
84 //
85 // sample the energy rate of the scattered gamma
86 //
87 G4double epsilon, epsilonsq, onecost, sint2, greject;
88
89 G4double epsilon0 = 1./(1. + 2*E0_m) , epsilon0sq = epsilon0*epsilon0;
90 G4double alpha1 = - log(epsilon0) , alpha2 = 0.5*(1.- epsilon0sq);
91
92 do {
93 if (alpha1/(alpha1+alpha2) > G4UniformRand())
94 { epsilon = exp(-alpha1*G4UniformRand()); // epsilon0**r
95 epsilonsq = epsilon*epsilon; }
96 else {
97 epsilonsq = epsilon0sq + (1.- epsilon0sq)*G4UniformRand();
98 epsilon = sqrt(epsilonsq);
99 };
100 onecost = (1.- epsilon)/(epsilon*E0_m);
101 sint2 = onecost*(2.-onecost);
102 greject = 1. - epsilon*sint2/(1.+ epsilonsq);
103 } while (greject < G4UniformRand());
104
105
106 //
107 // Phi determination
108 //
109 G4double minimum=0., maximum=twopi, middle=0., resolution=0.001;
110 G4double Rand = G4UniformRand();
111
112 int j = 0;
113 while ((j < 100) && (std::abs(SetPhi(epsilon,sint2,middle,Rand)) > resolution))
114 {
115 middle = (maximum + minimum)/2;
116 if (SetPhi(epsilon,sint2,middle,Rand)*
117 SetPhi(epsilon,sint2,minimum,Rand)<0) maximum = middle;
118 else minimum = middle;
119 j++;
120 }
121
122 //
123 // scattered gamma angles. ( Z - axis along the parent gamma)
124 //
125 G4double cosTeta = 1. - onecost , sinTeta = sqrt (sint2);
126 G4double Phi = middle;
127 G4double dirx = sinTeta*cos(Phi), diry = sinTeta*sin(Phi), dirz = cosTeta;
128
129 //
130 // update G4VParticleChange for the scattered gamma
131 //
132 G4double GammaEnergy1 = epsilon*GammaEnergy0;
133
134 // New polarization
135 //
136 G4ThreeVector GammaPolarization1 = SetNewPolarization(epsilon,sint2,Phi,
137 cosTeta,
138 GammaPolarization0);
139 // Set new direction
140 G4ThreeVector GammaDirection1 ( dirx,diry,dirz );
141
142 // Change reference frame.
143 SystemOfRefChange(GammaDirection0,GammaDirection1,
144 GammaPolarization0,GammaPolarization1);
145
146 G4double localEnergyDeposit = 0.;
147
148 if (GammaEnergy1 > fminimalEnergy)
149 {
150 aParticleChange.ProposeEnergy(GammaEnergy1);
151 }
152 else
153 {
154 localEnergyDeposit += GammaEnergy1;
155 aParticleChange.ProposeEnergy(0.) ;
156 aParticleChange.ProposeTrackStatus(fStopAndKill);
157 }
158
159 //
160 // kinematic of the scattered electron
161 //
162 G4double ElecKineEnergy = GammaEnergy0 - GammaEnergy1;
163
164 if (ElecKineEnergy > fminimalEnergy)
165 {
166 G4double ElecMomentum = sqrt(ElecKineEnergy*
167 (ElecKineEnergy+2.*electron_mass_c2));
168 G4ThreeVector ElecDirection (
169 (GammaEnergy0*GammaDirection0 - GammaEnergy1*GammaDirection1)
170 *(1./ElecMomentum));
171
172 // create G4DynamicParticle object for the electron.
173 G4DynamicParticle* aElectron= new G4DynamicParticle (
174 G4Electron::Electron(),ElecDirection,ElecKineEnergy);
175 aParticleChange.SetNumberOfSecondaries(1);
176 aParticleChange.AddSecondary( aElectron );
177 }
178 else
179 {
180 aParticleChange.SetNumberOfSecondaries(0);
181 localEnergyDeposit += ElecKineEnergy;
182 }
183
184 aParticleChange.ProposeLocalEnergyDeposit(localEnergyDeposit);
185
186 // Reset NbOfInteractionLengthLeft and return aParticleChange
187 return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep);
188}
189
190//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191
192G4double G4PolarizedComptonScattering::SetPhi(G4double EnergyRate,
193 G4double sinsqrth,
194 G4double phi,
195 G4double rand)
196{
197 G4double cosphi = cos(phi), sinphi = sin(phi);
198 G4double PhiDetermination = ((twopi*rand - phi)
199 *(EnergyRate + 1./EnergyRate - sinsqrth))
200 + (sinsqrth*sinphi*cosphi);
201 return PhiDetermination;
202}
203
204//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
205
206G4ThreeVector G4PolarizedComptonScattering::SetNewPolarization(
207 G4double EnergyRate,
208 G4double sinsqrth,
209 G4double phi,
210 G4double costheta,
211 G4ThreeVector&)
212{
213 G4double cosphi = cos(phi), sinphi = sin(phi);
214 //// G4double ParallelIntensityPolar = EnergyRate + 1./EnergyRate
215 //// + 2. - 4.*sinsqrth*cosphi*cosphi;
216 G4double ParallelIntensityPolar = EnergyRate + 1./EnergyRate
217 - 2.*sinsqrth*cosphi*cosphi;
218 G4double PerpendiIntensityPolar = EnergyRate + 1./EnergyRate - 2.;
219 G4double PolarizationDegree = sqrt(sinsqrth*sinphi*sinphi+costheta*costheta);
220 G4double sintheta = sqrt(sinsqrth);
221
222 G4ThreeVector GammaPolarization1;
223 // depolarization probability (1-P)
224 if ( G4UniformRand() > (PerpendiIntensityPolar/ParallelIntensityPolar) )
225 {
226 // Parallel to initial polarization
227 GammaPolarization1.setX(PolarizationDegree);
228 GammaPolarization1.setY(-sinsqrth*sinphi*cosphi/PolarizationDegree);
229 GammaPolarization1.setZ(-sintheta*costheta*cosphi/PolarizationDegree);
230
231 }
232 else
233 {
234 // Perpendicular to initial polarization
235 GammaPolarization1.setX(0.);
236 GammaPolarization1.setY(costheta/PolarizationDegree);
237 GammaPolarization1.setZ(-sintheta*sinphi/PolarizationDegree);
238
239 };
240
241 return GammaPolarization1;
242}
243
244//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
245
246void G4PolarizedComptonScattering::SystemOfRefChange(G4ThreeVector& Direction0,
247 G4ThreeVector& Direction1,
248 G4ThreeVector& Polarization0,
249 G4ThreeVector& Polarization1)
250{
251 // Angles for go back to the original RS
252 G4double cosTeta0 = Direction0.cosTheta(), sinTeta0 = sin(Direction0.theta());
253 G4double cosPhi0 = cos(Direction0.phi()), sinPhi0 = sin(Direction0.phi());
254
255
256
257 G4double cosPsi, sinPsi;
258
259 if (sinTeta0 != 0. )
260 {
261 cosPsi = -Polarization0.z()/sinTeta0;
262 if (cosPhi0 != 0.)
263 sinPsi = (Polarization0.y() - cosTeta0*sinPhi0*cosPsi)/cosPhi0;
264 else sinPsi = -Polarization0.x()/sinPhi0;
265
266 }
267 else
268 {
269 cosPsi = Polarization0.x()/cosTeta0;
270 sinPsi = Polarization0.y();
271 }
272 G4double Psi = atan(sinPsi/cosPsi);
273
274 // Rotation along Z axe
275 Direction1.rotateZ(Psi);
276 //
277 Direction1.rotateUz(Direction0);
278 aParticleChange.ProposeMomentumDirection(Direction1);
279
280 // 3 Euler angles rotation for scattered photon polarization
281 Polarization1.rotateZ(Psi);
282 Polarization1.rotateUz(Direction0);
283 aParticleChange.ProposePolarization(Polarization1);
284
285}
286
287//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.