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

Last change on this file since 1006 was 1006, checked in by garnier, 15 years ago

fichiers oublies
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-ref-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.