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