source: trunk/source/processes/electromagnetic/polarisation/src/G4PolarizedCompton.cc @ 1337

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

tag geant4.9.4 beta 1 + modifs locales
File size: 13.3 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.9 2008/10/30 22:34:23 schaelic Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-01 $
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  SetProcessSubType(fComptonScattering);
84}
85
86//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
87 
88G4PolarizedCompton::~G4PolarizedCompton()
89{
90  if (theAsymmetryTable) {
91    theAsymmetryTable->clearAndDestroy();
92    delete theAsymmetryTable;
93  }
94}
95
96//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97
98void G4PolarizedCompton::InitialiseProcess(const G4ParticleDefinition*)
99{
100  if(!isInitialised) {
101    isInitialised = true;
102    SetBuildTableFlag(true);
103    SetSecondaryParticle(G4Electron::Electron());
104    G4double emin = MinKinEnergy();
105    G4double emax = MaxKinEnergy();
106    emModel = new G4PolarizedComptonModel();
107    if(0 == mType) selectedModel = new G4KleinNishinaCompton();
108    else if(10 == mType) selectedModel = emModel;
109    selectedModel->SetLowEnergyLimit(emin);
110    selectedModel->SetHighEnergyLimit(emax);
111    AddEmModel(1, selectedModel);
112  } 
113}
114
115//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
116
117void G4PolarizedCompton::PrintInfo()
118{
119  G4cout << " Total cross sections has a good parametrisation"
120         << " from 10 KeV to (100/Z) GeV" 
121         << "\n      Sampling according " << selectedModel->GetName() << " model" 
122         << G4endl;
123}         
124
125//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
126
127void G4PolarizedCompton::SetModel(const G4String& s)
128{
129  if(s == "Klein-Nishina") mType = 0;
130  if(s == "Polarized-Compton") mType = 10;
131}
132
133
134//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
135//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
136
137
138
139G4double G4PolarizedCompton::GetMeanFreePath(
140                                   const G4Track& aTrack,
141                                   G4double   previousStepSize,
142                                   G4ForceCondition* condition)
143{
144  // *** get unploarised mean free path from lambda table ***
145  G4double mfp = G4VEmProcess::GetMeanFreePath(aTrack, previousStepSize, condition);
146
147
148   if (theAsymmetryTable && useAsymmetryTable) {
149     // *** get asymmetry, if target is polarized ***
150     const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
151     const G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
152     const G4StokesVector GammaPolarization = aTrack.GetPolarization();
153     const G4ParticleMomentum GammaDirection0 = aDynamicGamma->GetMomentumDirection();
154
155     G4Material*         aMaterial = aTrack.GetMaterial();
156     G4VPhysicalVolume*  aPVolume  = aTrack.GetVolume();
157     G4LogicalVolume*    aLVolume  = aPVolume->GetLogicalVolume();
158
159     //   G4Material* bMaterial = aLVolume->GetMaterial();
160     G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance();
161
162     const G4bool VolumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
163     G4StokesVector ElectronPolarization = polarizationManger->GetVolumePolarization(aLVolume);
164
165     if (!VolumeIsPolarized || mfp == DBL_MAX) return mfp;
166     
167     if (verboseLevel>=2) {
168
169       G4cout << " Mom " << GammaDirection0  << G4endl;
170       G4cout << " Polarization " << GammaPolarization  << G4endl;
171       G4cout << " MaterialPol. " << ElectronPolarization  << G4endl;
172       G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
173       G4cout << " Log. Volume  " << aLVolume->GetName() << G4endl;
174       G4cout << " Material     " << aMaterial          << G4endl;
175     }
176
177     G4int midx= CurrentMaterialCutsCoupleIndex();
178     G4PhysicsVector * aVector=(*theAsymmetryTable)(midx);
179     
180     G4double asymmetry=0;
181     if (aVector) {
182       G4bool isOutRange;
183       asymmetry = aVector->GetValue(GammaEnergy, isOutRange);
184     } else {
185       G4cout << " MaterialIndex     " << midx << " is out of range \n";
186       asymmetry=0;
187     }
188
189     //  we have to determine angle between particle motion
190     //  and target polarisation here 
191     //      circ pol * Vec(ElectronPol)*Vec(PhotonMomentum)
192     //  both vectors in global reference frame
193     
194     G4double pol=ElectronPolarization*GammaDirection0;
195     
196     G4double polProduct = GammaPolarization.p3() * pol;
197     mfp *= 1. / ( 1. + polProduct * asymmetry );
198
199     if (verboseLevel>=2) {
200       G4cout << " MeanFreePath:  " << mfp / mm << " mm " << G4endl;
201       G4cout << " Asymmetry:     " << asymmetry          << G4endl;
202       G4cout << " PolProduct:    " << polProduct         << G4endl;
203     }
204   }
205
206   return mfp;
207}
208
209G4double G4PolarizedCompton::PostStepGetPhysicalInteractionLength(
210                                   const G4Track& aTrack,
211                                   G4double   previousStepSize,
212                                   G4ForceCondition* condition)
213{
214  // *** get unploarised mean free path from lambda table ***
215  G4double mfp = G4VEmProcess::PostStepGetPhysicalInteractionLength(aTrack, previousStepSize, condition);
216
217
218   if (theAsymmetryTable && useAsymmetryTable) {
219     // *** get asymmetry, if target is polarized ***
220     const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
221     const G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
222     const G4StokesVector GammaPolarization = aTrack.GetPolarization();
223     const G4ParticleMomentum GammaDirection0 = aDynamicGamma->GetMomentumDirection();
224
225     G4Material*         aMaterial = aTrack.GetMaterial();
226     G4VPhysicalVolume*  aPVolume  = aTrack.GetVolume();
227     G4LogicalVolume*    aLVolume  = aPVolume->GetLogicalVolume();
228
229     //   G4Material* bMaterial = aLVolume->GetMaterial();
230     G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance();
231
232     const G4bool VolumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
233     G4StokesVector ElectronPolarization = polarizationManger->GetVolumePolarization(aLVolume);
234
235     if (!VolumeIsPolarized || mfp == DBL_MAX) return mfp;
236     
237     if (verboseLevel>=2) {
238
239       G4cout << " Mom " << GammaDirection0  << G4endl;
240       G4cout << " Polarization " << GammaPolarization  << G4endl;
241       G4cout << " MaterialPol. " << ElectronPolarization  << G4endl;
242       G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
243       G4cout << " Log. Volume  " << aLVolume->GetName() << G4endl;
244       G4cout << " Material     " << aMaterial          << G4endl;
245     }
246
247     G4int midx= CurrentMaterialCutsCoupleIndex();
248     G4PhysicsVector * aVector=(*theAsymmetryTable)(midx);
249     
250     G4double asymmetry=0;
251     if (aVector) {
252       G4bool isOutRange;
253       asymmetry = aVector->GetValue(GammaEnergy, isOutRange);
254     } else {
255       G4cout << " MaterialIndex     " << midx << " is out of range \n";
256       asymmetry=0;
257     }
258
259     //  we have to determine angle between particle motion
260     //  and target polarisation here 
261     //      circ pol * Vec(ElectronPol)*Vec(PhotonMomentum)
262     //  both vectors in global reference frame
263     
264     G4double pol=ElectronPolarization*GammaDirection0;
265     
266     G4double polProduct = GammaPolarization.p3() * pol;
267     mfp *= 1. / ( 1. + polProduct * asymmetry );
268
269     if (verboseLevel>=2) {
270       G4cout << " MeanFreePath:  " << mfp / mm << " mm " << G4endl;
271       G4cout << " Asymmetry:     " << asymmetry          << G4endl;
272       G4cout << " PolProduct:    " << polProduct         << G4endl;
273     }
274   }
275
276   return mfp;
277}
278
279//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
280
281void G4PolarizedCompton::PreparePhysicsTable(const G4ParticleDefinition& part)
282{
283  G4VEmProcess::PreparePhysicsTable(part);
284  if(buildAsymmetryTable)
285    theAsymmetryTable = G4PhysicsTableHelper::PreparePhysicsTable(theAsymmetryTable);
286}
287
288//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
289
290
291void G4PolarizedCompton::BuildPhysicsTable(const G4ParticleDefinition& part)
292{
293  // *** build (unpolarized) cross section tables (Lambda)
294  G4VEmProcess::BuildPhysicsTable(part);
295  if(buildAsymmetryTable)
296    BuildAsymmetryTable(part);
297}
298
299//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
300
301
302void G4PolarizedCompton::BuildAsymmetryTable(const G4ParticleDefinition& part)
303{
304  // Access to materials
305  const G4ProductionCutsTable* theCoupleTable=
306        G4ProductionCutsTable::GetProductionCutsTable();
307  size_t numOfCouples = theCoupleTable->GetTableSize();
308  for(size_t i=0; i<numOfCouples; ++i) {
309    if (!theAsymmetryTable) break;
310    if (theAsymmetryTable->GetFlag(i)) {
311
312      // create physics vector and fill it
313      const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
314      // use same parameters as for lambda
315      G4PhysicsVector* aVector = LambdaPhysicsVector(couple);
316      //      modelManager->FillLambdaVector(aVector, couple, startFromNull);
317
318      for (G4int j = 0 ; j < LambdaBinning() ; ++j ) {
319        G4double lowEdgeEnergy = aVector->GetLowEdgeEnergy(j);
320        G4double tasm=0.;
321        G4double asym = ComputeAsymmetry(lowEdgeEnergy, couple, part, 0., tasm);
322        aVector->PutValue(j,asym);
323      }
324
325      G4PhysicsTableHelper::SetPhysicsVector(theAsymmetryTable, i, aVector);
326    }
327  }
328
329}
330
331//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
332
333
334G4double G4PolarizedCompton::ComputeAsymmetry(G4double energy,
335                         const G4MaterialCutsCouple* couple,
336                       const G4ParticleDefinition& aParticle,
337                               G4double cut,
338                               G4double & tAsymmetry)
339{
340  G4double lAsymmetry = 0.0;
341  tAsymmetry=0;
342
343  //
344  // calculate polarized cross section
345  //
346  G4ThreeVector thePolarization=G4ThreeVector(0.,0.,1.);
347  emModel->SetTargetPolarization(thePolarization);
348  emModel->SetBeamPolarization(thePolarization);
349  G4double sigma2=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
350
351  //
352  // calculate unpolarized cross section
353  //
354  thePolarization=G4ThreeVector();
355  emModel->SetTargetPolarization(thePolarization);
356  emModel->SetBeamPolarization(thePolarization);
357  G4double sigma0=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
358
359  // determine assymmetries
360  if (sigma0>0.) {
361    lAsymmetry=sigma2/sigma0-1.;
362  }
363  return lAsymmetry;
364}
365
366//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Note: See TracBrowser for help on using the repository browser.