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

Last change on this file since 1350 was 1337, checked in by garnier, 14 years ago

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