source: trunk/source/processes/electromagnetic/polarisation/src/G4PolarizedComptonModel.cc@ 819

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

import all except CVS
File size: 13.7 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: G4PolarizedComptonModel.cc,v 1.4 2007/05/23 08:52:20 vnivanch Exp $
28// GEANT4 tag $Name: $
29//
30// -------------------------------------------------------------------
31//
32// GEANT4 Class file
33//
34//
35// File name: G4PolarizedComptonModel
36//
37// Author: Andreas Schaelicke
38//
39// Creation date: 01.05.2005
40//
41// Modifications:
42// 18-07-06 use newly calculated cross sections (P. Starovoitov)
43// 21-08-05 update interface (A. Schaelicke)
44//
45// Class Description:
46//
47// -------------------------------------------------------------------
48//
49//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51
52#include "G4PolarizedComptonModel.hh"
53#include "G4Electron.hh"
54#include "G4Gamma.hh"
55#include "Randomize.hh"
56#include "G4DataVector.hh"
57#include "G4ParticleChangeForGamma.hh"
58
59
60#include "G4StokesVector.hh"
61#include "G4PolarizationManager.hh"
62#include "G4PolarizationHelper.hh"
63#include "G4PolarizedComptonCrossSection.hh"
64
65//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66
67G4PolarizedComptonModel::G4PolarizedComptonModel(const G4ParticleDefinition*,
68 const G4String& nam)
69 : G4KleinNishinaCompton(0,nam),
70 verboseLevel(0)
71{
72 crossSectionCalculator=new G4PolarizedComptonCrossSection();
73}
74
75//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76
77G4PolarizedComptonModel::~G4PolarizedComptonModel()
78{
79 if (crossSectionCalculator) delete crossSectionCalculator;
80}
81
82
83
84G4double G4PolarizedComptonModel::ComputeAsymmetryPerAtom
85 (G4double gammaEnergy, G4double /*Z*/)
86
87{
88 G4double asymmetry = 0.0 ;
89
90 G4double k0 = gammaEnergy / electron_mass_c2 ;
91 G4double k1 = 1 + 2*k0 ;
92
93 asymmetry = -k0;
94 asymmetry *= (k0 + 1.)*sqr(k1)*std::log(k1) - 2.*k0*(5.*sqr(k0) + 4.*k0 + 1.);
95 asymmetry /= ((k0 - 2.)*k0 -2.)*sqr(k1)*std::log(k1) + 2.*k0*(k0*(k0 + 1.)*(k0 + 8.) + 2.);
96
97 // G4cout<<"energy = "<<GammaEnergy<<" asymmetry = "<<asymmetry<<"\t\t GAM = "<<k0<<G4endl;
98 if (asymmetry>1.) G4cout<<"ERROR in G4PolarizedComptonModel::ComputeAsymmetryPerAtom"<<G4endl;
99
100 return asymmetry;
101}
102
103
104G4double G4PolarizedComptonModel::ComputeCrossSectionPerAtom(
105 const G4ParticleDefinition* pd,
106 G4double kinEnergy,
107 G4double Z,
108 G4double A,
109 G4double cut,
110 G4double emax)
111{
112 double xs =
113 G4KleinNishinaCompton::ComputeCrossSectionPerAtom(pd,kinEnergy,
114 Z,A,cut,emax);
115 G4double polzz = theBeamPolarization.p3()*theTargetPolarization.z();
116 if (polzz!=0) {
117 G4double asym=ComputeAsymmetryPerAtom(kinEnergy, Z);
118 xs*=(1.+polzz*asym);
119 }
120 return xs;
121}
122
123
124//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
125
126void G4PolarizedComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
127 const G4MaterialCutsCouple*,
128 const G4DynamicParticle* aDynamicGamma,
129 G4double,
130 G4double)
131{
132 const G4Track * aTrack = fParticleChange->GetCurrentTrack();
133 G4VPhysicalVolume* aPVolume = aTrack->GetVolume();
134 G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
135
136 if (verboseLevel>=1)
137 G4cout<<"G4PolarizedComptonModel::SampleSecondaries in "
138 << aLVolume->GetName() <<G4endl;
139
140 G4PolarizationManager * polarizationManager = G4PolarizationManager::GetInstance();
141
142 // obtain polarization of the beam
143 theBeamPolarization = aDynamicGamma->GetPolarization();
144 theBeamPolarization.SetPhoton();
145
146 // obtain polarization of the media
147 const G4bool targetIsPolarized = polarizationManager->IsPolarized(aLVolume);
148 theTargetPolarization = polarizationManager->GetVolumePolarization(aLVolume);
149
150 // if beam is linear polarized or target is transversely polarized
151 // determine the angle to x-axis
152 // (assumes same PRF as in the polarization definition)
153
154 G4ThreeVector gamDirection0 = aDynamicGamma->GetMomentumDirection();
155
156 // transfere theTargetPolarization
157 // into the gamma frame (problem electron is at rest)
158 if (targetIsPolarized)
159 theTargetPolarization.rotateUz(gamDirection0);
160
161 // The scattered gamma energy is sampled according to Klein - Nishina formula.
162 // The random number techniques of Butcher & Messel are used
163 // (Nuc Phys 20(1960),15).
164 // Note : Effects due to binding of atomic electrons are negliged.
165
166 G4double gamEnergy0 = aDynamicGamma->GetKineticEnergy();
167 G4double E0_m = gamEnergy0 / electron_mass_c2 ;
168
169
170 //
171 // sample the energy rate of the scattered gamma
172 //
173
174 G4double epsilon, epsilonsq, onecost, sint2, greject ;
175
176 G4double epsilon0 = 1./(1. + 2.*E0_m);
177 G4double epsilon0sq = epsilon0*epsilon0;
178 G4double alpha1 = - std::log(epsilon0);
179 G4double alpha2 = 0.5*(1.- epsilon0sq);
180
181 G4double polarization = theBeamPolarization.p3()*theTargetPolarization.p3();
182 do {
183 if ( alpha1/(alpha1+alpha2) > G4UniformRand() ) {
184 epsilon = std::exp(-alpha1*G4UniformRand()); // epsilon0**r
185 epsilonsq = epsilon*epsilon;
186
187 } else {
188 epsilonsq = epsilon0sq + (1.- epsilon0sq)*G4UniformRand();
189 epsilon = std::sqrt(epsilonsq);
190 };
191
192 onecost = (1.- epsilon)/(epsilon*E0_m);
193 sint2 = onecost*(2.-onecost);
194
195
196 G4double gdiced = 2.*(1./epsilon+epsilon);
197 G4double gdist = 1./epsilon + epsilon - sint2
198 - polarization*(1./epsilon-epsilon)*(1.-onecost);
199
200 greject = gdist/gdiced;
201
202 if (greject>1) G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n"
203 <<" costh rejection does not work properly: "<<greject<<G4endl;
204
205 } while (greject < G4UniformRand());
206
207 //
208 // scattered gamma angles. ( Z - axis along the parent gamma)
209 //
210
211 G4double cosTeta = 1. - onecost;
212 G4double sinTeta = std::sqrt (sint2);
213 G4double Phi;
214 do {
215 Phi = twopi * G4UniformRand();
216 G4double gdiced = 1./epsilon + epsilon - sint2
217 + std::abs(theBeamPolarization.p3())*
218 ( std::abs((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3())
219 +(1.-epsilon)*sinTeta*(std::sqrt(sqr(theTargetPolarization.p1())
220 + sqr(theTargetPolarization.p2()))))
221 +sint2*(std::sqrt(sqr(theBeamPolarization.p1()) + sqr(theBeamPolarization.p2())));
222
223 G4double gdist = 1./epsilon + epsilon - sint2
224 + theBeamPolarization.p3()*
225 ((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3()
226 +(1.-epsilon)*sinTeta*(std::cos(Phi)*theTargetPolarization.p1()+
227 std::sin(Phi)*theTargetPolarization.p2()))
228 -sint2*(std::cos(2.*Phi)*theBeamPolarization.p1()
229 +std::sin(2.*Phi)*theBeamPolarization.p2());
230 greject = gdist/gdiced;
231
232 if (greject>1.+1.e-10 || greject<0) G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n"
233 <<" phi rejection does not work properly: "<<greject<<G4endl;
234
235 if (greject<1.e-3) {
236 G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n"
237 <<" phi rejection does not work properly: "<<greject<<"\n";
238 G4cout<<" greject="<<greject<<" phi="<<Phi<<" cost="<<cosTeta<<"\n";
239 G4cout<<" gdiced="<<gdiced<<" gdist="<<gdist<<"\n";
240 G4cout<<" eps="<<epsilon<<" 1/eps="<<1./epsilon<<"\n";
241 }
242
243 } while (greject < G4UniformRand());
244 G4double dirx = sinTeta*std::cos(Phi), diry = sinTeta*std::sin(Phi), dirz = cosTeta;
245
246 //
247 // update G4VParticleChange for the scattered gamma
248 //
249
250 G4ThreeVector gamDirection1 ( dirx,diry,dirz );
251 gamDirection1.rotateUz(gamDirection0);
252 G4double gamEnergy1 = epsilon*gamEnergy0;
253 fParticleChange->SetProposedKineticEnergy(gamEnergy1);
254
255
256 if(gamEnergy1 > lowestGammaEnergy) {
257 fParticleChange->ProposeMomentumDirection(gamDirection1);
258 } else {
259 fParticleChange->ProposeTrackStatus(fStopAndKill);
260 gamEnergy1 += fParticleChange->GetLocalEnergyDeposit();
261 fParticleChange->ProposeLocalEnergyDeposit(gamEnergy1);
262 }
263
264 //
265 // kinematic of the scattered electron
266 //
267
268 G4double eKinEnergy = gamEnergy0 - gamEnergy1;
269 G4ThreeVector eDirection = gamEnergy0*gamDirection0 - gamEnergy1*gamDirection1;
270 eDirection = eDirection.unit();
271
272 //
273 // calculate Stokesvector of final state photon and electron
274 //
275 G4ThreeVector nInteractionFrame;
276 if((gamEnergy1 > lowestGammaEnergy) ||
277 (eKinEnergy > DBL_MIN)) {
278
279 // determine interaction plane
280// nInteractionFrame =
281// G4PolarizationHelper::GetFrame(gamDirection1,eDirection);
282 if (gamEnergy1 > lowestGammaEnergy)
283 nInteractionFrame = G4PolarizationHelper::GetFrame(gamDirection1,gamDirection0);
284 else
285 nInteractionFrame = G4PolarizationHelper::GetFrame(gamDirection0, eDirection);
286
287 // transfere theBeamPolarization and theTargetPolarization
288 // into the interaction frame (note electron is in gamma frame)
289 if (verboseLevel>=1) {
290 G4cout << "========================================\n";
291 G4cout << " nInteractionFrame = " <<nInteractionFrame<<"\n";
292 G4cout << " GammaDirection0 = " <<gamDirection0<<"\n";
293 G4cout << " gammaPolarization = " <<theBeamPolarization<<"\n";
294 G4cout << " electronPolarization = " <<theTargetPolarization<<"\n";
295 }
296
297 theBeamPolarization.InvRotateAz(nInteractionFrame,gamDirection0);
298 theTargetPolarization.InvRotateAz(nInteractionFrame,gamDirection0);
299
300 if (verboseLevel>=1) {
301 G4cout << "----------------------------------------\n";
302 G4cout << " gammaPolarization = " <<theBeamPolarization<<"\n";
303 G4cout << " electronPolarization = " <<theTargetPolarization<<"\n";
304 G4cout << "----------------------------------------\n";
305 }
306
307 // initialize the polarization transfer matrix
308 crossSectionCalculator->Initialize(epsilon,E0_m,0.,
309 theBeamPolarization,
310 theTargetPolarization,2);
311 }
312
313 // if(eKinEnergy > DBL_MIN)
314 {
315 // in interaction frame
316 // calculate polarization transfer to the photon (in interaction plane)
317 finalGammaPolarization = crossSectionCalculator->GetPol2();
318 if (verboseLevel>=1) G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n";
319 finalGammaPolarization.SetPhoton();
320
321 // translate polarization into particle reference frame
322 finalGammaPolarization.RotateAz(nInteractionFrame,gamDirection1);
323 //store polarization vector
324 fParticleChange->ProposePolarization(finalGammaPolarization);
325 if (finalGammaPolarization.mag() > 1.+1.e-8){
326 G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl;
327 G4cout<<"Polarization of final photon more than 100%"<<G4endl;
328 G4cout<<finalGammaPolarization<<" mag = "<<finalGammaPolarization.mag()<<G4endl;
329 }
330 if (verboseLevel>=1) {
331 G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n";
332 G4cout << " GammaDirection1 = " <<gamDirection1<<"\n";
333 }
334 }
335
336 // if (ElecKineEnergy > fminimalEnergy) {
337 {
338 finalElectronPolarization = crossSectionCalculator->GetPol3();
339 if (verboseLevel>=1)
340 G4cout << " electronPolarization1 = " <<finalElectronPolarization<<"\n";
341
342 // transfer into particle reference frame
343 finalElectronPolarization.RotateAz(nInteractionFrame,eDirection);
344 if (verboseLevel>=1) {
345 G4cout << " electronPolarization1 = " <<finalElectronPolarization<<"\n";
346 G4cout << " ElecDirection = " <<eDirection<<"\n";
347 }
348 }
349 if (verboseLevel>=1)
350 G4cout << "========================================\n";
351
352
353 if(eKinEnergy > DBL_MIN) {
354
355 // create G4DynamicParticle object for the electron.
356 G4DynamicParticle* aElectron = new G4DynamicParticle(theElectron,eDirection,eKinEnergy);
357 //store polarization vector
358 if (finalElectronPolarization.mag() > 1.+1.e-8){
359 G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl;
360 G4cout<<"Polarization of final electron more than 100%"<<G4endl;
361 G4cout<<finalElectronPolarization<<" mag = "<<finalElectronPolarization.mag()<<G4endl;
362 }
363 aElectron->SetPolarization(finalElectronPolarization.p1(),
364 finalElectronPolarization.p2(),
365 finalElectronPolarization.p3());
366 fvect->push_back(aElectron);
367 }
368}
369
370//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
371
372
Note: See TracBrowser for help on using the repository browser.