source: trunk/source/processes/electromagnetic/polarisation/src/G4PolarizedCompton.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.2 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: G4PolarizedCompton.cc,v 1.7 2007/07/10 09:35:37 schaelic Exp $
28// GEANT4 tag $Name: geant4-09-01-patch-02 $
29//
30//
31// File name: G4PolarizedCompton
32//
33// Author: Andreas Schaelicke
34// based on code by Michel Maire / Vladimir IVANTCHENKO
35// Class description
36//
37// modified version respecting media and beam polarization
38// using the stokes formalism
39//
40// Creation date: 01.05.2005
41//
42// Modifications:
43//
44// 01-01-05, include polarization description (A.Stahl)
45// 01-01-05, create asymmetry table and determine interactionlength (A.Stahl)
46// 01-05-05, update handling of media polarization (A.Schalicke)
47// 01-05-05, update polarized differential cross section (A.Schalicke)
48// 20-05-05, added polarization transfer (A.Schalicke)
49// 10-06-05, transformation between different reference frames (A.Schalicke)
50// 17-10-05, correct reference frame dependence in GetMeanFreePath (A.Schalicke)
51// 26-07-06, cross section recalculated (P.Starovoitov)
52// 09-08-06, make it work under current geant4 release (A.Schalicke)
53// 11-06-07, add PostStepGetPhysicalInteractionLength (A.Schalicke)
54// -----------------------------------------------------------------------------
55
56
57#include "G4PolarizedCompton.hh"
58#include "G4Electron.hh"
59
60#include "G4StokesVector.hh"
61#include "G4PolarizationManager.hh"
62#include "G4PolarizedComptonModel.hh"
63#include "G4ProductionCutsTable.hh"
64#include "G4PhysicsTableHelper.hh"
65#include "G4KleinNishinaCompton.hh"
66#include "G4PolarizedComptonModel.hh"
67
68//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69
70G4PolarizedCompton::G4PolarizedCompton(const G4String& processName,
71 G4ProcessType type):
72 G4VEmProcess (processName, type),
73 buildAsymmetryTable(true),
74 useAsymmetryTable(true),
75 isInitialised(false),
76 selectedModel(0),
77 mType(10),
78 theAsymmetryTable(NULL)
79{
80 SetLambdaBinning(90);
81 SetMinKinEnergy(0.1*keV);
82 SetMaxKinEnergy(100.0*GeV);
83}
84
85//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
86
87G4PolarizedCompton::~G4PolarizedCompton()
88{
89 if (theAsymmetryTable) {
90 theAsymmetryTable->clearAndDestroy();
91 delete theAsymmetryTable;
92 }
93}
94
95//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96
97void G4PolarizedCompton::InitialiseProcess(const G4ParticleDefinition*)
98{
99 if(!isInitialised) {
100 isInitialised = true;
101 SetBuildTableFlag(true);
102 SetSecondaryParticle(G4Electron::Electron());
103 G4double emin = MinKinEnergy();
104 G4double emax = MaxKinEnergy();
105 emModel = new G4PolarizedComptonModel();
106 if(0 == mType) selectedModel = new G4KleinNishinaCompton();
107 else if(10 == mType) selectedModel = emModel;
108 selectedModel->SetLowEnergyLimit(emin);
109 selectedModel->SetHighEnergyLimit(emax);
110 AddEmModel(1, selectedModel);
111 }
112}
113
114//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
115
116void G4PolarizedCompton::PrintInfo()
117{
118 G4cout << " Total cross sections has a good parametrisation"
119 << " from 10 KeV to (100/Z) GeV"
120 << "\n Sampling according " << selectedModel->GetName() << " model"
121 << G4endl;
122}
123
124//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
125
126void G4PolarizedCompton::SetModel(const G4String& s)
127{
128 if(s == "Klein-Nishina") mType = 0;
129 if(s == "Polarized-Compton") mType = 10;
130}
131
132
133//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
134//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
135
136
137
138G4double G4PolarizedCompton::GetMeanFreePath(
139 const G4Track& aTrack,
140 G4double previousStepSize,
141 G4ForceCondition* condition)
142{
143 // *** get unploarised mean free path from lambda table ***
144 G4double mfp = G4VEmProcess::GetMeanFreePath(aTrack, previousStepSize, condition);
145
146
147 if (theAsymmetryTable && useAsymmetryTable) {
148 // *** get asymmetry, if target is polarized ***
149 const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
150 const G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
151 const G4StokesVector GammaPolarization = aTrack.GetPolarization();
152 const G4ParticleMomentum GammaDirection0 = aDynamicGamma->GetMomentumDirection();
153
154 G4Material* aMaterial = aTrack.GetMaterial();
155 G4VPhysicalVolume* aPVolume = aTrack.GetVolume();
156 G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
157
158 // G4Material* bMaterial = aLVolume->GetMaterial();
159 G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance();
160
161 const G4bool VolumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
162 G4StokesVector ElectronPolarization = polarizationManger->GetVolumePolarization(aLVolume);
163
164 if (!VolumeIsPolarized || mfp == DBL_MAX) return mfp;
165
166 if (verboseLevel>=2) {
167
168 G4cout << " Mom " << GammaDirection0 << G4endl;
169 G4cout << " Polarization " << GammaPolarization << G4endl;
170 G4cout << " MaterialPol. " << ElectronPolarization << G4endl;
171 G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
172 G4cout << " Log. Volume " << aLVolume->GetName() << G4endl;
173 G4cout << " Material " << aMaterial << G4endl;
174 }
175
176 G4int midx= CurrentMaterialCutsCoupleIndex();
177 G4PhysicsVector * aVector=(*theAsymmetryTable)(midx);
178
179 G4double asymmetry=0;
180 if (aVector) {
181 G4bool isOutRange;
182 asymmetry = aVector->GetValue(GammaEnergy, isOutRange);
183 } else {
184 G4cout << " MaterialIndex " << midx << " is out of range \n";
185 asymmetry=0;
186 }
187
188 // we have to determine angle between particle motion
189 // and target polarisation here
190 // circ pol * Vec(ElectronPol)*Vec(PhotonMomentum)
191 // both vectors in global reference frame
192
193 G4double pol=ElectronPolarization*GammaDirection0;
194
195 G4double polProduct = GammaPolarization.p3() * pol;
196 mfp *= 1. / ( 1. + polProduct * asymmetry );
197
198 if (verboseLevel>=2) {
199 G4cout << " MeanFreePath: " << mfp / mm << " mm " << G4endl;
200 G4cout << " Asymmetry: " << asymmetry << G4endl;
201 G4cout << " PolProduct: " << polProduct << G4endl;
202 }
203 }
204
205 return mfp;
206}
207
208G4double G4PolarizedCompton::PostStepGetPhysicalInteractionLength(
209 const G4Track& aTrack,
210 G4double previousStepSize,
211 G4ForceCondition* condition)
212{
213 // *** get unploarised mean free path from lambda table ***
214 G4double mfp = G4VEmProcess::PostStepGetPhysicalInteractionLength(aTrack, previousStepSize, condition);
215
216
217 if (theAsymmetryTable && useAsymmetryTable) {
218 // *** get asymmetry, if target is polarized ***
219 const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
220 const G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
221 const G4StokesVector GammaPolarization = aTrack.GetPolarization();
222 const G4ParticleMomentum GammaDirection0 = aDynamicGamma->GetMomentumDirection();
223
224 G4Material* aMaterial = aTrack.GetMaterial();
225 G4VPhysicalVolume* aPVolume = aTrack.GetVolume();
226 G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
227
228 // G4Material* bMaterial = aLVolume->GetMaterial();
229 G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance();
230
231 const G4bool VolumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
232 G4StokesVector ElectronPolarization = polarizationManger->GetVolumePolarization(aLVolume);
233
234 if (!VolumeIsPolarized || mfp == DBL_MAX) return mfp;
235
236 if (verboseLevel>=2) {
237
238 G4cout << " Mom " << GammaDirection0 << G4endl;
239 G4cout << " Polarization " << GammaPolarization << G4endl;
240 G4cout << " MaterialPol. " << ElectronPolarization << G4endl;
241 G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
242 G4cout << " Log. Volume " << aLVolume->GetName() << G4endl;
243 G4cout << " Material " << aMaterial << G4endl;
244 }
245
246 G4int midx= CurrentMaterialCutsCoupleIndex();
247 G4PhysicsVector * aVector=(*theAsymmetryTable)(midx);
248
249 G4double asymmetry=0;
250 if (aVector) {
251 G4bool isOutRange;
252 asymmetry = aVector->GetValue(GammaEnergy, isOutRange);
253 } else {
254 G4cout << " MaterialIndex " << midx << " is out of range \n";
255 asymmetry=0;
256 }
257
258 // we have to determine angle between particle motion
259 // and target polarisation here
260 // circ pol * Vec(ElectronPol)*Vec(PhotonMomentum)
261 // both vectors in global reference frame
262
263 G4double pol=ElectronPolarization*GammaDirection0;
264
265 G4double polProduct = GammaPolarization.p3() * pol;
266 mfp *= 1. / ( 1. + polProduct * asymmetry );
267
268 if (verboseLevel>=2) {
269 G4cout << " MeanFreePath: " << mfp / mm << " mm " << G4endl;
270 G4cout << " Asymmetry: " << asymmetry << G4endl;
271 G4cout << " PolProduct: " << polProduct << G4endl;
272 }
273 }
274
275 return mfp;
276}
277
278//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
279
280void G4PolarizedCompton::PreparePhysicsTable(const G4ParticleDefinition& part)
281{
282 G4VEmProcess::PreparePhysicsTable(part);
283 if(buildAsymmetryTable)
284 theAsymmetryTable = G4PhysicsTableHelper::PreparePhysicsTable(theAsymmetryTable);
285}
286
287//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
288
289
290void G4PolarizedCompton::BuildPhysicsTable(const G4ParticleDefinition& part)
291{
292 // *** build (unpolarized) cross section tables (Lambda)
293 G4VEmProcess::BuildPhysicsTable(part);
294 if(buildAsymmetryTable)
295 BuildAsymmetryTable(part);
296}
297
298//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
299
300
301void G4PolarizedCompton::BuildAsymmetryTable(const G4ParticleDefinition& part)
302{
303 // Access to materials
304 const G4ProductionCutsTable* theCoupleTable=
305 G4ProductionCutsTable::GetProductionCutsTable();
306 size_t numOfCouples = theCoupleTable->GetTableSize();
307 for(size_t i=0; i<numOfCouples; ++i) {
308 if (!theAsymmetryTable) break;
309 if (theAsymmetryTable->GetFlag(i)) {
310
311 // create physics vector and fill it
312 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
313 // use same parameters as for lambda
314 G4PhysicsVector* aVector = LambdaPhysicsVector(couple);
315 // modelManager->FillLambdaVector(aVector, couple, startFromNull);
316
317 for (G4int j = 0 ; j < LambdaBinning() ; ++j ) {
318 G4double lowEdgeEnergy = aVector->GetLowEdgeEnergy(j);
319 G4double tasm=0.;
320 G4double asym = ComputeAsymmetry(lowEdgeEnergy, couple, part, 0., tasm);
321 aVector->PutValue(j,asym);
322 }
323
324 G4PhysicsTableHelper::SetPhysicsVector(theAsymmetryTable, i, aVector);
325 }
326 }
327
328}
329
330//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
331
332
333G4double G4PolarizedCompton::ComputeAsymmetry(G4double energy,
334 const G4MaterialCutsCouple* couple,
335 const G4ParticleDefinition& aParticle,
336 G4double cut,
337 G4double & tAsymmetry)
338{
339 G4double lAsymmetry = 0.0;
340 tAsymmetry=0;
341
342 //
343 // calculate polarized cross section
344 //
345 G4ThreeVector thePolarization=G4ThreeVector(0.,0.,1.);
346 emModel->SetTargetPolarization(thePolarization);
347 emModel->SetBeamPolarization(thePolarization);
348 G4double sigma2=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
349
350 //
351 // calculate unpolarized cross section
352 //
353 thePolarization=G4ThreeVector();
354 emModel->SetTargetPolarization(thePolarization);
355 emModel->SetBeamPolarization(thePolarization);
356 G4double sigma0=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
357
358 // determine assymmetries
359 if (sigma0>0.) {
360 lAsymmetry=sigma2/sigma0-1.;
361 }
362 return lAsymmetry;
363}
364
365//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Note: See TracBrowser for help on using the repository browser.