Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PolarizedCompton.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 //
27 // $Id: G4PolarizedCompton.cc 69847 2013-05-16 09:36:18Z gcosmo $
28 //
29 //
30 // File name: G4PolarizedCompton
31 //
32 // Author: Andreas Schaelicke
33 // based on code by Michel Maire / Vladimir IVANTCHENKO
34 // Class description
35 //
36 // modified version respecting media and beam polarization
37 // using the stokes formalism
38 //
39 // Creation date: 01.05.2005
40 //
41 // Modifications:
42 //
43 // 01-01-05, include polarization description (A.Stahl)
44 // 01-01-05, create asymmetry table and determine interactionlength (A.Stahl)
45 // 01-05-05, update handling of media polarization (A.Schalicke)
46 // 01-05-05, update polarized differential cross section (A.Schalicke)
47 // 20-05-05, added polarization transfer (A.Schalicke)
48 // 10-06-05, transformation between different reference frames (A.Schalicke)
49 // 17-10-05, correct reference frame dependence in GetMeanFreePath (A.Schalicke)
50 // 26-07-06, cross section recalculated (P.Starovoitov)
51 // 09-08-06, make it work under current geant4 release (A.Schalicke)
52 // 11-06-07, add PostStepGetPhysicalInteractionLength (A.Schalicke)
53 // -----------------------------------------------------------------------------
54 
55 
56 #include "G4PolarizedCompton.hh"
57 #include "G4SystemOfUnits.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 
70 G4PolarizedCompton::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  emModel = 0;
85 }
86 
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88 
89 G4PolarizedCompton::~G4PolarizedCompton()
90 {
91  if (theAsymmetryTable) {
92  theAsymmetryTable->clearAndDestroy();
93  delete theAsymmetryTable;
94  }
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
98 
99 void G4PolarizedCompton::InitialiseProcess(const G4ParticleDefinition*)
100 {
101  if(!isInitialised) {
102  isInitialised = true;
103  SetBuildTableFlag(true);
104  SetSecondaryParticle(G4Electron::Electron());
105  G4double emin = MinKinEnergy();
106  G4double emax = MaxKinEnergy();
107  emModel = new G4PolarizedComptonModel();
108  if(0 == mType) selectedModel = new G4KleinNishinaCompton();
109  else if(10 == mType) selectedModel = emModel;
110  selectedModel->SetLowEnergyLimit(emin);
111  selectedModel->SetHighEnergyLimit(emax);
112  AddEmModel(1, selectedModel);
113  }
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117 
118 void G4PolarizedCompton::PrintInfo()
119 {
120  G4cout << " Total cross sections has a good parametrisation"
121  << " from 10 KeV to (100/Z) GeV"
122  << "\n Sampling according " << selectedModel->GetName() << " model"
123  << G4endl;
124 }
125 
126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
127 
128 void G4PolarizedCompton::SetModel(const G4String& ss)
129 {
130  if(ss == "Klein-Nishina") mType = 0;
131  if(ss == "Polarized-Compton") mType = 10;
132 }
133 
134 
135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
137 
138 
139 
140 G4double G4PolarizedCompton::GetMeanFreePath(
141  const G4Track& aTrack,
142  G4double previousStepSize,
143  G4ForceCondition* condition)
144 {
145  // *** get unploarised mean free path from lambda table ***
146  G4double mfp = G4VEmProcess::GetMeanFreePath(aTrack, previousStepSize, condition);
147 
148 
149  if (theAsymmetryTable && useAsymmetryTable) {
150  // *** get asymmetry, if target is polarized ***
151  const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
152  const G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
153  const G4StokesVector GammaPolarization = aTrack.GetPolarization();
154  const G4ParticleMomentum GammaDirection0 = aDynamicGamma->GetMomentumDirection();
155 
156  G4Material* aMaterial = aTrack.GetMaterial();
157  G4VPhysicalVolume* aPVolume = aTrack.GetVolume();
158  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
159 
160  // G4Material* bMaterial = aLVolume->GetMaterial();
161  G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance();
162 
163  const G4bool VolumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
164  G4StokesVector ElectronPolarization = polarizationManger->GetVolumePolarization(aLVolume);
165 
166  if (!VolumeIsPolarized || mfp == DBL_MAX) return mfp;
167 
168  if (verboseLevel>=2) {
169 
170  G4cout << " Mom " << GammaDirection0 << G4endl;
171  G4cout << " Polarization " << GammaPolarization << G4endl;
172  G4cout << " MaterialPol. " << ElectronPolarization << G4endl;
173  G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
174  G4cout << " Log. Volume " << aLVolume->GetName() << G4endl;
175  G4cout << " Material " << aMaterial << G4endl;
176  }
177 
178  G4int midx= CurrentMaterialCutsCoupleIndex();
179  G4PhysicsVector * aVector=(*theAsymmetryTable)(midx);
180 
181  G4double asymmetry=0;
182  if (aVector) {
183  G4bool isOutRange;
184  asymmetry = aVector->GetValue(GammaEnergy, isOutRange);
185  } else {
186  G4cout << " MaterialIndex " << midx << " is out of range \n";
187  asymmetry=0;
188  }
189 
190  // we have to determine angle between particle motion
191  // and target polarisation here
192  // circ pol * Vec(ElectronPol)*Vec(PhotonMomentum)
193  // both vectors in global reference frame
194 
195  G4double pol=ElectronPolarization*GammaDirection0;
196 
197  G4double polProduct = GammaPolarization.p3() * pol;
198  mfp *= 1. / ( 1. + polProduct * asymmetry );
199 
200  if (verboseLevel>=2) {
201  G4cout << " MeanFreePath: " << mfp / mm << " mm " << G4endl;
202  G4cout << " Asymmetry: " << asymmetry << G4endl;
203  G4cout << " PolProduct: " << polProduct << G4endl;
204  }
205  }
206 
207  return mfp;
208 }
209 
210 G4double G4PolarizedCompton::PostStepGetPhysicalInteractionLength(
211  const G4Track& aTrack,
212  G4double previousStepSize,
213  G4ForceCondition* condition)
214 {
215  // *** get unploarised mean free path from lambda table ***
216  G4double mfp = G4VEmProcess::PostStepGetPhysicalInteractionLength(aTrack, previousStepSize, condition);
217 
218 
219  if (theAsymmetryTable && useAsymmetryTable) {
220  // *** get asymmetry, if target is polarized ***
221  const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
222  const G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
223  const G4StokesVector GammaPolarization = aTrack.GetPolarization();
224  const G4ParticleMomentum GammaDirection0 = aDynamicGamma->GetMomentumDirection();
225 
226  G4Material* aMaterial = aTrack.GetMaterial();
227  G4VPhysicalVolume* aPVolume = aTrack.GetVolume();
228  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
229 
230  // G4Material* bMaterial = aLVolume->GetMaterial();
231  G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance();
232 
233  const G4bool VolumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
234  G4StokesVector ElectronPolarization = polarizationManger->GetVolumePolarization(aLVolume);
235 
236  if (!VolumeIsPolarized || mfp == DBL_MAX) return mfp;
237 
238  if (verboseLevel>=2) {
239 
240  G4cout << " Mom " << GammaDirection0 << G4endl;
241  G4cout << " Polarization " << GammaPolarization << G4endl;
242  G4cout << " MaterialPol. " << ElectronPolarization << G4endl;
243  G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
244  G4cout << " Log. Volume " << aLVolume->GetName() << G4endl;
245  G4cout << " Material " << aMaterial << G4endl;
246  }
247 
248  G4int midx= CurrentMaterialCutsCoupleIndex();
249  G4PhysicsVector * aVector=(*theAsymmetryTable)(midx);
250 
251  G4double asymmetry=0;
252  if (aVector) {
253  G4bool isOutRange;
254  asymmetry = aVector->GetValue(GammaEnergy, isOutRange);
255  } else {
256  G4cout << " MaterialIndex " << midx << " is out of range \n";
257  asymmetry=0;
258  }
259 
260  // we have to determine angle between particle motion
261  // and target polarisation here
262  // circ pol * Vec(ElectronPol)*Vec(PhotonMomentum)
263  // both vectors in global reference frame
264 
265  G4double pol=ElectronPolarization*GammaDirection0;
266 
267  G4double polProduct = GammaPolarization.p3() * pol;
268  mfp *= 1. / ( 1. + polProduct * asymmetry );
269 
270  if (verboseLevel>=2) {
271  G4cout << " MeanFreePath: " << mfp / mm << " mm " << G4endl;
272  G4cout << " Asymmetry: " << asymmetry << G4endl;
273  G4cout << " PolProduct: " << polProduct << G4endl;
274  }
275  }
276 
277  return mfp;
278 }
279 
280 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
281 
282 void G4PolarizedCompton::PreparePhysicsTable(const G4ParticleDefinition& part)
283 {
284  G4VEmProcess::PreparePhysicsTable(part);
285  if(buildAsymmetryTable)
286  theAsymmetryTable = G4PhysicsTableHelper::PreparePhysicsTable(theAsymmetryTable);
287 }
288 
289 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
290 
291 
292 void G4PolarizedCompton::BuildPhysicsTable(const G4ParticleDefinition& part)
293 {
294  // *** build (unpolarized) cross section tables (Lambda)
295  G4VEmProcess::BuildPhysicsTable(part);
296  if(buildAsymmetryTable)
297  BuildAsymmetryTable(part);
298 }
299 
300 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
301 
302 
303 void G4PolarizedCompton::BuildAsymmetryTable(const G4ParticleDefinition& part)
304 {
305  // Access to materials
306  const G4ProductionCutsTable* theCoupleTable=
307  G4ProductionCutsTable::GetProductionCutsTable();
308  size_t numOfCouples = theCoupleTable->GetTableSize();
309  for(size_t i=0; i<numOfCouples; ++i) {
310  if (!theAsymmetryTable) break;
311  if (theAsymmetryTable->GetFlag(i)) {
312 
313  // create physics vector and fill it
314  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
315  // use same parameters as for lambda
316  G4PhysicsVector* aVector = LambdaPhysicsVector(couple);
317  // modelManager->FillLambdaVector(aVector, couple, startFromNull);
318 
319  for (G4int j = 0 ; j < LambdaBinning() ; ++j ) {
320  G4double lowEdgeEnergy = aVector->GetLowEdgeEnergy(j);
321  G4double tasm=0.;
322  G4double asym = ComputeAsymmetry(lowEdgeEnergy, couple, part, 0., tasm);
323  aVector->PutValue(j,asym);
324  }
325 
326  G4PhysicsTableHelper::SetPhysicsVector(theAsymmetryTable, i, aVector);
327  }
328  }
329 
330 }
331 
332 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
333 
334 
335 G4double G4PolarizedCompton::ComputeAsymmetry(G4double energy,
336  const G4MaterialCutsCouple* couple,
337  const G4ParticleDefinition& aParticle,
338  G4double cut,
339  G4double & tAsymmetry)
340 {
341  G4double lAsymmetry = 0.0;
342  tAsymmetry=0;
343 
344  //
345  // calculate polarized cross section
346  //
347  G4ThreeVector thePolarization=G4ThreeVector(0.,0.,1.);
348  emModel->SetTargetPolarization(thePolarization);
349  emModel->SetBeamPolarization(thePolarization);
350  G4double sigma2=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
351 
352  //
353  // calculate unpolarized cross section
354  //
355  thePolarization=G4ThreeVector();
356  emModel->SetTargetPolarization(thePolarization);
357  emModel->SetBeamPolarization(thePolarization);
358  G4double sigma0=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
359 
360  // determine assymmetries
361  if (sigma0>0.) {
362  lAsymmetry=sigma2/sigma0-1.;
363  }
364  return lAsymmetry;
365 }
366 
367 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Generated on Sat May 25 2013 14:33:36 for Geant4 by   doxygen 1.8.4