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

Last change on this file since 900 was 819, checked in by garnier, 17 years ago

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