source: trunk/source/processes/electromagnetic/polarisation/src/G4eplusPolarizedAnnihilation.cc@ 900

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

import all except CVS
File size: 14.7 KB
RevLine 
[819]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// $Id: G4eplusPolarizedAnnihilation.cc,v 1.6 2007/10/02 11:36:44 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-01-patch-02 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4eplusPolarizedAnnihilation
35//
36// Author: A. Schaelicke on base of Vladimir Ivanchenko / Michel Maire code
37//
38// Creation date: 02.07.2006
39//
40// Modifications:
41// 26-07-06 modified cross section (P. Starovoitov)
42// 21-08-06 interface updated (A. Schaelicke)
43// 11-06-07, add PostStepGetPhysicalInteractionLength (A.Schalicke)
44// 02-10-07, enable AtRest (V.Ivanchenko)
45//
46//
47// Class Description:
48//
49// Polarized process of e+ annihilation into 2 gammas
50//
51
52//
53// -------------------------------------------------------------------
54//
55//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
57
58#include "G4eplusPolarizedAnnihilation.hh"
59#include "G4MaterialCutsCouple.hh"
60#include "G4Gamma.hh"
61#include "G4PhysicsVector.hh"
62#include "G4PhysicsLogVector.hh"
63
64
65#include "G4PolarizedAnnihilationModel.hh"
66#include "G4PhysicsTableHelper.hh"
67#include "G4ProductionCutsTable.hh"
68#include "G4PolarizationManager.hh"
69#include "G4PolarizationHelper.hh"
70#include "G4StokesVector.hh"
71
72//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73
74G4eplusPolarizedAnnihilation::G4eplusPolarizedAnnihilation(const G4String& name)
75 : G4VEmProcess(name), isInitialised(false),
76 theAsymmetryTable(NULL),
77 theTransverseAsymmetryTable(NULL)
78{
79 enableAtRestDoIt = true;
80}
81
82//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83
84G4eplusPolarizedAnnihilation::~G4eplusPolarizedAnnihilation()
85{
86 if (theAsymmetryTable) {
87 theAsymmetryTable->clearAndDestroy();
88 delete theAsymmetryTable;
89 }
90 if (theTransverseAsymmetryTable) {
91 theTransverseAsymmetryTable->clearAndDestroy();
92 delete theTransverseAsymmetryTable;
93 }
94}
95
96//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97
98void G4eplusPolarizedAnnihilation::InitialiseProcess(const G4ParticleDefinition*)
99{
100 if(!isInitialised) {
101 isInitialised = true;
102 // SetVerboseLevel(3);
103 SetBuildTableFlag(true);
104 SetStartFromNullFlag(false);
105 SetSecondaryParticle(G4Gamma::Gamma());
106 G4double emin = 0.1*keV;
107 G4double emax = 100.*TeV;
108 SetLambdaBinning(120);
109 SetMinKinEnergy(emin);
110 SetMaxKinEnergy(emax);
111 emModel = new G4PolarizedAnnihilationModel();
112 emModel->SetLowEnergyLimit(emin);
113 emModel->SetHighEnergyLimit(emax);
114 AddEmModel(1, emModel);
115 }
116}
117
118
119//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
120
121 // for polarization
122
123G4double G4eplusPolarizedAnnihilation::GetMeanFreePath(const G4Track& track,
124 G4double previousStepSize,
125 G4ForceCondition* condition)
126{
127 G4double mfp = G4VEmProcess::GetMeanFreePath(track, previousStepSize, condition);
128
129 if (theAsymmetryTable) {
130
131 G4Material* aMaterial = track.GetMaterial();
132 G4VPhysicalVolume* aPVolume = track.GetVolume();
133 G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
134
135 // G4Material* bMaterial = aLVolume->GetMaterial();
136 G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance();
137
138 const G4bool volumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
139 G4StokesVector electronPolarization = polarizationManger->GetVolumePolarization(aLVolume);
140
141 if (!volumeIsPolarized || mfp == DBL_MAX) return mfp;
142
143 // *** get asymmetry, if target is polarized ***
144 const G4DynamicParticle* aDynamicPositron = track.GetDynamicParticle();
145 const G4double positronEnergy = aDynamicPositron->GetKineticEnergy();
146 const G4StokesVector positronPolarization = track.GetPolarization();
147 const G4ParticleMomentum positronDirection0 = aDynamicPositron->GetMomentumDirection();
148
149 if (verboseLevel>=2) {
150
151 G4cout << " Mom " << positronDirection0 << G4endl;
152 G4cout << " Polarization " << positronPolarization << G4endl;
153 G4cout << " MaterialPol. " << electronPolarization << G4endl;
154 G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
155 G4cout << " Log. Volume " << aLVolume->GetName() << G4endl;
156 G4cout << " Material " << aMaterial << G4endl;
157 }
158
159 G4bool isOutRange;
160 G4int idx= CurrentMaterialCutsCoupleIndex();
161 G4double lAsymmetry = (*theAsymmetryTable)(idx)->
162 GetValue(positronEnergy, isOutRange);
163 G4double tAsymmetry = (*theTransverseAsymmetryTable)(idx)->
164 GetValue(positronEnergy, isOutRange);
165
166 G4double polZZ = positronPolarization.z()*
167 electronPolarization*positronDirection0;
168 G4double polXX = positronPolarization.x()*
169 electronPolarization*G4PolarizationHelper::GetParticleFrameX(positronDirection0);
170 G4double polYY = positronPolarization.y()*
171 electronPolarization*G4PolarizationHelper::GetParticleFrameY(positronDirection0);
172
173 G4double impact = 1. + polZZ*lAsymmetry + (polXX + polYY)*tAsymmetry;
174
175 mfp *= 1. / impact;
176
177 if (verboseLevel>=2) {
178 G4cout << " MeanFreePath: " << mfp / mm << " mm " << G4endl;
179 G4cout << " Asymmetry: " << lAsymmetry << ", " << tAsymmetry << G4endl;
180 G4cout << " PolProduct: " << polXX << ", " << polYY << ", " << polZZ << G4endl;
181 }
182 }
183
184 return mfp;
185}
186
187//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
188
189G4double G4eplusPolarizedAnnihilation::PostStepGetPhysicalInteractionLength(
190 const G4Track& track,
191 G4double previousStepSize,
192 G4ForceCondition* condition)
193{
194 G4double mfp = G4VEmProcess::PostStepGetPhysicalInteractionLength(track, previousStepSize, condition);
195
196 if (theAsymmetryTable) {
197
198 G4Material* aMaterial = track.GetMaterial();
199 G4VPhysicalVolume* aPVolume = track.GetVolume();
200 G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
201
202 // G4Material* bMaterial = aLVolume->GetMaterial();
203 G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance();
204
205 const G4bool volumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
206 G4StokesVector electronPolarization = polarizationManger->GetVolumePolarization(aLVolume);
207
208 if (!volumeIsPolarized || mfp == DBL_MAX) return mfp;
209
210 // *** get asymmetry, if target is polarized ***
211 const G4DynamicParticle* aDynamicPositron = track.GetDynamicParticle();
212 const G4double positronEnergy = aDynamicPositron->GetKineticEnergy();
213 const G4StokesVector positronPolarization = track.GetPolarization();
214 const G4ParticleMomentum positronDirection0 = aDynamicPositron->GetMomentumDirection();
215
216 if (verboseLevel>=2) {
217
218 G4cout << " Mom " << positronDirection0 << G4endl;
219 G4cout << " Polarization " << positronPolarization << G4endl;
220 G4cout << " MaterialPol. " << electronPolarization << G4endl;
221 G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
222 G4cout << " Log. Volume " << aLVolume->GetName() << G4endl;
223 G4cout << " Material " << aMaterial << G4endl;
224 }
225
226 G4bool isOutRange;
227 G4int idx= CurrentMaterialCutsCoupleIndex();
228 G4double lAsymmetry = (*theAsymmetryTable)(idx)->
229 GetValue(positronEnergy, isOutRange);
230 G4double tAsymmetry = (*theTransverseAsymmetryTable)(idx)->
231 GetValue(positronEnergy, isOutRange);
232
233 G4double polZZ = positronPolarization.z()*
234 electronPolarization*positronDirection0;
235 G4double polXX = positronPolarization.x()*
236 electronPolarization*G4PolarizationHelper::GetParticleFrameX(positronDirection0);
237 G4double polYY = positronPolarization.y()*
238 electronPolarization*G4PolarizationHelper::GetParticleFrameY(positronDirection0);
239
240 G4double impact = 1. + polZZ*lAsymmetry + (polXX + polYY)*tAsymmetry;
241
242 mfp *= 1. / impact;
243
244 if (verboseLevel>=2) {
245 G4cout << " MeanFreePath: " << mfp / mm << " mm " << G4endl;
246 G4cout << " Asymmetry: " << lAsymmetry << ", " << tAsymmetry << G4endl;
247 G4cout << " PolProduct: " << polXX << ", " << polYY << ", " << polZZ << G4endl;
248 }
249 }
250
251 return mfp;
252}
253//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
254
255void G4eplusPolarizedAnnihilation::BuildPhysicsTable(const G4ParticleDefinition& pd)
256{
257 G4VEmProcess::BuildPhysicsTable(pd);
258 BuildAsymmetryTable(pd);
259}
260//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
261
262void G4eplusPolarizedAnnihilation::PreparePhysicsTable(const G4ParticleDefinition& pd)
263{
264 G4VEmProcess::PreparePhysicsTable(pd);
265 theAsymmetryTable = G4PhysicsTableHelper::PreparePhysicsTable(theAsymmetryTable);
266 theTransverseAsymmetryTable = G4PhysicsTableHelper::PreparePhysicsTable(theTransverseAsymmetryTable);
267}
268
269void G4eplusPolarizedAnnihilation::BuildAsymmetryTable(const G4ParticleDefinition& part)
270{
271 // Access to materials
272 const G4ProductionCutsTable* theCoupleTable=
273 G4ProductionCutsTable::GetProductionCutsTable();
274 size_t numOfCouples = theCoupleTable->GetTableSize();
275 G4cout<<" annih-numOfCouples="<<numOfCouples<<"\n";
276 for(size_t i=0; i<numOfCouples; ++i) {
277 G4cout<<"annih- "<<i<<"/"<<numOfCouples<<"\n";
278 if (!theAsymmetryTable) break;
279 G4cout<<"annih- "<<theAsymmetryTable->GetFlag(i)<<"\n";
280 if (theAsymmetryTable->GetFlag(i)) {
281 G4cout<<" building pol-annih ... \n";
282
283 // create physics vector and fill it
284 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
285
286 // use same parameters as for lambda
287 G4PhysicsVector* aVector = LambdaPhysicsVector(couple);
288 G4PhysicsVector* tVector = LambdaPhysicsVector(couple);
289
290 for (G4int j = 0 ; j < LambdaBinning() ; ++j ) {
291 G4double lowEdgeEnergy = aVector->GetLowEdgeEnergy(j);
292 G4double tasm=0.;
293 G4double asym = ComputeAsymmetry(lowEdgeEnergy, couple, part, 0., tasm);
294 aVector->PutValue(j,asym);
295 tVector->PutValue(j,tasm);
296 }
297
298 G4PhysicsTableHelper::SetPhysicsVector(theAsymmetryTable, i, aVector);
299 G4PhysicsTableHelper::SetPhysicsVector(theTransverseAsymmetryTable, i, tVector);
300 }
301 }
302
303}
304
305
306//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
307
308G4double G4eplusPolarizedAnnihilation::ComputeAsymmetry(G4double energy,
309 const G4MaterialCutsCouple* couple,
310 const G4ParticleDefinition& aParticle,
311 G4double cut,
312 G4double &tAsymmetry)
313{
314 G4double lAsymmetry = 0.0;
315 tAsymmetry = 0.0;
316
317 // calculate polarized cross section
318 theTargetPolarization=G4ThreeVector(0.,0.,1.);
319 emModel->SetTargetPolarization(theTargetPolarization);
320 emModel->SetBeamPolarization(theTargetPolarization);
321 G4double sigma2=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
322
323 // calculate transversely polarized cross section
324 theTargetPolarization=G4ThreeVector(1.,0.,0.);
325 emModel->SetTargetPolarization(theTargetPolarization);
326 emModel->SetBeamPolarization(theTargetPolarization);
327 G4double sigma3=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
328
329 // calculate unpolarized cross section
330 theTargetPolarization=G4ThreeVector();
331 emModel->SetTargetPolarization(theTargetPolarization);
332 emModel->SetBeamPolarization(theTargetPolarization);
333 G4double sigma0=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
334
335 // determine assymmetries
336 if (sigma0>0.) {
337 lAsymmetry=sigma2/sigma0-1.;
338 tAsymmetry=sigma3/sigma0-1.;
339 }
340 return lAsymmetry;
341
342}
343
344
345
346//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
347
348void G4eplusPolarizedAnnihilation::PrintInfo()
349{
350 G4cout << " Polarized model for annihilation into 2 photons"
351 << G4endl;
352}
353
354//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
355
356G4VParticleChange* G4eplusPolarizedAnnihilation::AtRestDoIt(const G4Track& aTrack,
357 const G4Step& )
358//
359// Performs the e+ e- annihilation when both particles are assumed at rest.
360// It generates two back to back photons with energy = electron_mass.
361// The angular distribution is isotropic.
362// GEANT4 internal units
363//
364// Note : Effects due to binding of atomic electrons are negliged.
365{
366 fParticleChange.InitializeForPostStep(aTrack);
367
368 fParticleChange.SetNumberOfSecondaries(2);
369
370 G4double cosTeta = 2.*G4UniformRand()-1. , sinTeta = std::sqrt(1.-cosTeta*cosTeta);
371 G4double phi = twopi * G4UniformRand();
372 G4ThreeVector direction (sinTeta*std::cos(phi), sinTeta*std::sin(phi), cosTeta);
373 fParticleChange.AddSecondary( new G4DynamicParticle (G4Gamma::Gamma(),
374 direction, electron_mass_c2) );
375 fParticleChange.AddSecondary( new G4DynamicParticle (G4Gamma::Gamma(),
376 -direction, electron_mass_c2) );
377 // Kill the incident positron
378 //
379 fParticleChange.ProposeTrackStatus(fStopAndKill);
380 return &fParticleChange;
381}
382
383//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Note: See TracBrowser for help on using the repository browser.