source: trunk/source/processes/electromagnetic/polarisation/src/G4ePolarizedIonisation.cc@ 880

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

import all except CVS
File size: 12.9 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: G4ePolarizedIonisation.cc,v 1.5 2007/07/10 09:35:37 schaelic Exp $
28// GEANT4 tag $Name: geant4-09-01-patch-02 $
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4ePolarizedIonisation
35//
36// Author: A.Schaelicke on base of Vladimir Ivanchenko code
37//
38// Creation date: 10.11.2005
39//
40// Modifications:
41//
42// 10-11-05, include polarization description (A.Schaelicke)
43// , create asymmetry table and determine interactionlength
44// , update polarized differential cross section
45//
46// 20-08-06, modified interface (A.Schaelicke)
47// 11-06-07, add PostStepGetPhysicalInteractionLength (A.Schalicke)
48//
49// Class Description:
50//
51//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
52//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53
54#include "G4ePolarizedIonisation.hh"
55#include "G4Electron.hh"
56#include "G4UniversalFluctuation.hh"
57#include "G4BohrFluctuations.hh"
58#include "G4UnitsTable.hh"
59
60#include "G4PolarizedMollerBhabhaModel.hh"
61#include "G4ProductionCutsTable.hh"
62#include "G4PolarizationManager.hh"
63#include "G4PolarizationHelper.hh"
64#include "G4StokesVector.hh"
65
66//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67
68G4ePolarizedIonisation::G4ePolarizedIonisation(const G4String& name)
69 : G4VEnergyLossProcess(name),
70 theElectron(G4Electron::Electron()),
71 isElectron(true),
72 isInitialised(false),
73 theAsymmetryTable(NULL),
74 theTransverseAsymmetryTable(NULL)
75{
76 verboseLevel=0;
77 SetDEDXBinning(120);
78 SetLambdaBinning(120);
79 numBinAsymmetryTable=120;
80
81 SetMinKinEnergy(0.1*keV);
82 SetMaxKinEnergy(100.0*TeV);
83 // PrintInfoDefinition();
84}
85
86//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
87
88G4ePolarizedIonisation::~G4ePolarizedIonisation()
89{
90 if (theAsymmetryTable) {
91 theAsymmetryTable->clearAndDestroy();
92 delete theAsymmetryTable;
93 }
94 if (theTransverseAsymmetryTable) {
95 theTransverseAsymmetryTable->clearAndDestroy();
96 delete theTransverseAsymmetryTable;
97 }
98}
99
100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101
102void G4ePolarizedIonisation::InitialiseEnergyLossProcess(
103 const G4ParticleDefinition* part,
104 const G4ParticleDefinition* /*part2*/)
105{
106 if(!isInitialised) {
107
108 if(part == G4Positron::Positron()) isElectron = false;
109 SetSecondaryParticle(theElectron);
110
111
112
113 flucModel = new G4UniversalFluctuation();
114 //flucModel = new G4BohrFluctuations();
115
116 // G4VEmModel* em = new G4MollerBhabhaModel();
117 emModel = new G4PolarizedMollerBhabhaModel;
118 emModel->SetLowEnergyLimit(100*eV);
119 emModel->SetHighEnergyLimit(100*TeV);
120 AddEmModel(1, emModel, flucModel);
121
122 isInitialised = true;
123 }
124}
125
126//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
127
128void G4ePolarizedIonisation::PrintInfo()
129{
130 G4cout << " Delta cross sections from Moller+Bhabha, "
131 << "good description from 1 KeV to 100 GeV."
132 << G4endl;
133}
134
135//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
136
137G4double G4ePolarizedIonisation::GetMeanFreePath(const G4Track& track,
138 G4double s,
139 G4ForceCondition* cond)
140{
141 // *** get unploarised mean free path from lambda table ***
142 G4double mfp = G4VEnergyLossProcess::GetMeanFreePath(track, s, cond);
143
144
145 // *** get asymmetry, if target is polarized ***
146 G4VPhysicalVolume* aPVolume = track.GetVolume();
147 G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
148
149 G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance();
150 const G4bool volumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
151 const G4StokesVector ePolarization = track.GetPolarization();
152
153 if (mfp != DBL_MAX && volumeIsPolarized && !ePolarization.IsZero()) {
154 const G4DynamicParticle* aDynamicElectron = track.GetDynamicParticle();
155 const G4double eEnergy = aDynamicElectron->GetKineticEnergy();
156 const G4ParticleMomentum eDirection0 = aDynamicElectron->GetMomentumDirection();
157
158 G4StokesVector volumePolarization = polarizationManger->GetVolumePolarization(aLVolume);
159
160 G4bool isOutRange;
161 size_t idx = CurrentMaterialCutsCoupleIndex();
162 G4double lAsymmetry = (*theAsymmetryTable)(idx)->
163 GetValue(eEnergy, isOutRange);
164 G4double tAsymmetry = (*theTransverseAsymmetryTable)(idx)->
165 GetValue(eEnergy, isOutRange);
166
167 // calculate longitudinal spin component
168 G4double polZZ = ePolarization.z()*
169 volumePolarization*eDirection0;
170 // calculate transvers spin components
171 G4double polXX = ePolarization.x()*
172 volumePolarization*G4PolarizationHelper::GetParticleFrameX(eDirection0);
173 G4double polYY = ePolarization.y()*
174 volumePolarization*G4PolarizationHelper::GetParticleFrameY(eDirection0);
175
176
177 G4double impact = 1. + polZZ*lAsymmetry + (polXX + polYY)*tAsymmetry;
178 // determine polarization dependent mean free path
179 mfp /= impact;
180 if (mfp <=0.) {
181 G4cout <<"PV impact ( "<<polXX<<" , "<<polYY<<" , "<<polZZ<<" )"<<G4endl;
182 G4cout << " impact on MFP is "<< impact <<G4endl;
183 G4cout<<" lAsymmetry= "<<lAsymmetry<<" ("<<std::fabs(lAsymmetry)-1.<<")\n";
184 G4cout<<" tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
185 }
186 }
187
188 return mfp;
189}
190
191G4double G4ePolarizedIonisation::PostStepGetPhysicalInteractionLength(const G4Track& track,
192 G4double s,
193 G4ForceCondition* cond)
194{
195 // *** get unploarised mean free path from lambda table ***
196 G4double mfp = G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength(track, s, cond);
197
198
199 // *** get asymmetry, if target is polarized ***
200 G4VPhysicalVolume* aPVolume = track.GetVolume();
201 G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
202
203 G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance();
204 const G4bool volumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
205 const G4StokesVector ePolarization = track.GetPolarization();
206
207 if (mfp != DBL_MAX && volumeIsPolarized && !ePolarization.IsZero()) {
208 const G4DynamicParticle* aDynamicElectron = track.GetDynamicParticle();
209 const G4double eEnergy = aDynamicElectron->GetKineticEnergy();
210 const G4ParticleMomentum eDirection0 = aDynamicElectron->GetMomentumDirection();
211
212 G4StokesVector volumePolarization = polarizationManger->GetVolumePolarization(aLVolume);
213
214 G4bool isOutRange;
215 size_t idx = CurrentMaterialCutsCoupleIndex();
216 G4double lAsymmetry = (*theAsymmetryTable)(idx)->
217 GetValue(eEnergy, isOutRange);
218 G4double tAsymmetry = (*theTransverseAsymmetryTable)(idx)->
219 GetValue(eEnergy, isOutRange);
220
221 // calculate longitudinal spin component
222 G4double polZZ = ePolarization.z()*
223 volumePolarization*eDirection0;
224 // calculate transvers spin components
225 G4double polXX = ePolarization.x()*
226 volumePolarization*G4PolarizationHelper::GetParticleFrameX(eDirection0);
227 G4double polYY = ePolarization.y()*
228 volumePolarization*G4PolarizationHelper::GetParticleFrameY(eDirection0);
229
230
231 G4double impact = 1. + polZZ*lAsymmetry + (polXX + polYY)*tAsymmetry;
232 // determine polarization dependent mean free path
233 mfp /= impact;
234 if (mfp <=0.) {
235 G4cout <<"PV impact ( "<<polXX<<" , "<<polYY<<" , "<<polZZ<<" )"<<G4endl;
236 G4cout << " impact on MFP is "<< impact <<G4endl;
237 G4cout<<" lAsymmetry= "<<lAsymmetry<<" ("<<std::fabs(lAsymmetry)-1.<<")\n";
238 G4cout<<" tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
239 }
240 }
241
242 return mfp;
243}
244
245//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
246void G4ePolarizedIonisation::BuildPhysicsTable(const G4ParticleDefinition& part)
247{
248 // *** build DEDX and (unpolarized) cross section tables
249 G4VEnergyLossProcess::BuildPhysicsTable(part);
250 // G4PhysicsTable* pt =
251 // BuildDEDXTable();
252
253
254 // *** build asymmetry-table
255 if (theAsymmetryTable) {
256 theAsymmetryTable->clearAndDestroy(); delete theAsymmetryTable;}
257 if (theTransverseAsymmetryTable) {
258 theTransverseAsymmetryTable->clearAndDestroy(); delete theTransverseAsymmetryTable;}
259
260 const G4ProductionCutsTable* theCoupleTable=
261 G4ProductionCutsTable::GetProductionCutsTable();
262 size_t numOfCouples = theCoupleTable->GetTableSize();
263
264 theAsymmetryTable = new G4PhysicsTable(numOfCouples);
265 theTransverseAsymmetryTable = new G4PhysicsTable(numOfCouples);
266
267 for (size_t j=0 ; j < numOfCouples; j++ ) {
268 // get cut value
269 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
270
271 G4double tcutmin = emModel->MinEnergyCut(&part, couple);
272 G4double cut = (*theCoupleTable->GetEnergyCutsVector(1))[j];
273 cut = std::max(cut, tcutmin);
274
275 //create physics vectors then fill it (same parameters as lambda vector)
276 G4PhysicsVector * ptrVectorA = LambdaPhysicsVector(couple,cut);
277 G4PhysicsVector * ptrVectorB = LambdaPhysicsVector(couple,cut);
278
279 for (size_t i = 0 ; i < numBinAsymmetryTable ; i++ ) {
280 G4double lowEdgeEnergy = ptrVectorA->GetLowEdgeEnergy(i);
281 G4double tasm=0.;
282 G4double asym = ComputeAsymmetry(lowEdgeEnergy, couple, part, cut, tasm);
283 ptrVectorA->PutValue(i,asym);
284 ptrVectorB->PutValue(i,tasm);
285 }
286 theAsymmetryTable->insertAt( j , ptrVectorA ) ;
287 theTransverseAsymmetryTable->insertAt( j , ptrVectorB ) ;
288 }
289
290}
291//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
292
293G4double G4ePolarizedIonisation::ComputeAsymmetry(G4double energy,
294 const G4MaterialCutsCouple* couple,
295 const G4ParticleDefinition& aParticle,
296 G4double cut,
297 G4double & tAsymmetry)
298{
299 G4double lAsymmetry = 0.0;
300 tAsymmetry = 0.0;
301 if (isElectron) {lAsymmetry = tAsymmetry = -1.0;}
302
303 // calculate polarized cross section
304 theTargetPolarization=G4ThreeVector(0.,0.,1.);
305 emModel->SetTargetPolarization(theTargetPolarization);
306 emModel->SetBeamPolarization(theTargetPolarization);
307 G4double sigma2=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
308
309 // calculate transversely polarized cross section
310 theTargetPolarization=G4ThreeVector(1.,0.,0.);
311 emModel->SetTargetPolarization(theTargetPolarization);
312 emModel->SetBeamPolarization(theTargetPolarization);
313 G4double sigma3=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
314
315 // calculate unpolarized cross section
316 theTargetPolarization=G4ThreeVector();
317 emModel->SetTargetPolarization(theTargetPolarization);
318 emModel->SetBeamPolarization(theTargetPolarization);
319 G4double sigma0=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
320 // determine assymmetries
321 if (sigma0>0.) {
322 lAsymmetry=sigma2/sigma0-1.;
323 tAsymmetry=sigma3/sigma0-1.;
324 }
325 if (std::fabs(lAsymmetry)>1.) {
326 G4cout<<" energy="<<energy<<"\n";
327 G4cout<<"WARNING lAsymmetry= "<<lAsymmetry<<" ("<<std::fabs(lAsymmetry)-1.<<")\n";
328 }
329 if (std::fabs(tAsymmetry)>1.) {
330 G4cout<<" energy="<<energy<<"\n";
331 G4cout<<"WARNING tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
332 }
333// else {
334// G4cout<<" tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
335// }
336 return lAsymmetry;
337}
338
339
Note: See TracBrowser for help on using the repository browser.