Geant4_10
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 76244 2013-11-08 11:12:59Z 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"
63 #include "G4ProductionCutsTable.hh"
64 #include "G4PhysicsTableHelper.hh"
65 #include "G4KleinNishinaCompton.hh"
67 
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69 
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);
84  emModel = 0;
85 }
86 
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88 
90 {
91  if (theAsymmetryTable) {
92  delete theAsymmetryTable;
93  }
94 }
95 
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97 
99 {
100  if(!isInitialised) {
101  isInitialised = true;
102  SetBuildTableFlag(true);
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 
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 
128 {
129  if(ss == "Klein-Nishina") mType = 0;
130  if(ss == "Polarized-Compton") mType = 10;
131 }
132 
133 
134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
136 
137 
138 
140  const G4Track& aTrack,
141  G4double previousStepSize,
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();
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 
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 
210  const G4Track& aTrack,
211  G4double previousStepSize,
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();
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 
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 
282 {
284  if(buildAsymmetryTable)
285  theAsymmetryTable = G4PhysicsTableHelper::PreparePhysicsTable(theAsymmetryTable);
286 }
287 
288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
289 
290 
292 {
293  // *** build (unpolarized) cross section tables (Lambda)
295  if(buildAsymmetryTable)
296  BuildAsymmetryTable(part);
297 }
298 
299 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
300 
301 
303 {
304  // Access to materials
305  const G4ProductionCutsTable* theCoupleTable=
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 
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....
G4double condition(const G4ErrorSymMatrix &m)
const G4ThreeVector & GetPolarization() const
static void SetPhysicsVector(G4PhysicsTable *physTable, size_t idx, G4PhysicsVector *vec)
G4double GetValue(G4double theEnergy, G4bool &isOutRange) const
G4int verboseLevel
Definition: G4VProcess.hh:368
G4String GetName() const
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
const G4DynamicParticle * GetDynamicParticle() const
void SetBuildTableFlag(G4bool val)
void BuildAsymmetryTable(const G4ParticleDefinition &part)
static G4PolarizationManager * GetInstance()
virtual void PreparePhysicsTable(const G4ParticleDefinition &)
G4double MaxKinEnergy() const
void BuildPhysicsTable(const G4ParticleDefinition &)
G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
G4double GetLowEdgeEnergy(size_t binNumber) const
int G4int
Definition: G4Types.hh:78
static G4PhysicsTable * PreparePhysicsTable(G4PhysicsTable *physTable)
G4int LambdaBinning() const
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
G4double p3() const
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:683
virtual G4double GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)
void SetLambdaBinning(G4int nbins)
size_t CurrentMaterialCutsCoupleIndex() const
G4GLOB_DLL std::ostream G4cout
double energy
Definition: plottest35.C:25
const G4String & GetName() const
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
G4double CrossSection(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:467
TString part[npart]
Definition: Style.C:32
void PutValue(size_t index, G4double theValue)
G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
void PreparePhysicsTable(const G4ParticleDefinition &)
G4Material * GetMaterial() const
static G4ProductionCutsTable * GetProductionCutsTable()
G4PolarizedCompton(const G4String &processName="pol-compt", G4ProcessType type=fElectromagnetic)
G4PhysicsVector * LambdaPhysicsVector(const G4MaterialCutsCouple *)
G4double ComputeAsymmetry(G4double energy, const G4MaterialCutsCouple *couple, const G4ParticleDefinition &particle, G4double cut, G4double &tAsymmetry)
void SetMaxKinEnergy(G4double e)
G4LogicalVolume * GetLogicalVolume() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
void SetBeamPolarization(const G4ThreeVector &pBeam)
void AddEmModel(G4int, G4VEmModel *, const G4Region *region=0)
void SetSecondaryParticle(const G4ParticleDefinition *p)
bool IsPolarized(G4LogicalVolume *lVol) const
G4VPhysicalVolume * GetVolume() const
virtual void InitialiseProcess(const G4ParticleDefinition *)
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
const G4String & GetName() const
Definition: G4VEmModel.hh:753
void SetMinKinEnergy(G4double e)
G4bool GetFlag(size_t i) const
double G4double
Definition: G4Types.hh:76
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:690
G4ForceCondition
G4double MinKinEnergy() const
#define DBL_MAX
Definition: templates.hh:83
const G4ThreeVector & GetVolumePolarization(G4LogicalVolume *lVol) const
void SetTargetPolarization(const G4ThreeVector &pTarget)
virtual void PrintInfo()
void SetModel(const G4String &name)
G4ProcessType