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

Last change on this file since 846 was 819, checked in by garnier, 16 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.