Geant4  10.01.p02
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 85018 2014-10-23 09:51:37Z 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 #include "G4EmParameters.hh"
68 
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
70 
72 
74  G4ProcessType type):
75  G4VEmProcess (processName, type),
76  buildAsymmetryTable(true),
77  useAsymmetryTable(true),
78  isInitialised(false),
79  mType(10)
80 {
82  SetBuildTableFlag(true);
86  SetSplineFlag(true);
87  emModel = 0;
88 }
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
91 
93 {
94  delete theAsymmetryTable;
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
98 
100 {
101  if(!isInitialised) {
102  isInitialised = true;
103  if(0 == mType) {
104  if(!EmModel(1)) { SetEmModel(new G4KleinNishinaCompton(), 1); }
105  } else {
107  SetEmModel(emModel, 1);
108  }
110  EmModel(1)->SetLowEnergyLimit(param->MinKinEnergy());
112  AddEmModel(1, EmModel(1));
113  }
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117 
119 {
120  G4cout << " Total cross sections has a good parametrisation"
121  << " from 10 KeV to (100/Z) GeV"
122  << "\n Sampling according " << EmModel(1)->GetName() << " model"
123  << G4endl;
124 }
125 
126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
127 
129 {
130  if(ss == "Klein-Nishina") { mType = 0; }
131  if(ss == "Polarized-Compton") { mType = 10; }
132 }
133 
134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
135 
137  const G4Track& aTrack,
138  G4double previousStepSize,
140 {
141  // *** get unploarised mean free path from lambda table ***
142  G4double mfp = G4VEmProcess::GetMeanFreePath(aTrack, previousStepSize, condition);
143 
145  // *** get asymmetry, if target is polarized ***
146  const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
147  const G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
148  const G4StokesVector GammaPolarization = aTrack.GetPolarization();
149  const G4ParticleMomentum GammaDirection0 = aDynamicGamma->GetMomentumDirection();
150 
151  G4Material* aMaterial = aTrack.GetMaterial();
152  G4VPhysicalVolume* aPVolume = aTrack.GetVolume();
153  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
154 
155  // G4Material* bMaterial = aLVolume->GetMaterial();
157 
158  G4bool VolumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
159  G4StokesVector ElectronPolarization = polarizationManger->GetVolumePolarization(aLVolume);
160 
161  if (!VolumeIsPolarized || mfp == DBL_MAX) return mfp;
162 
163  if (verboseLevel>=2) {
164 
165  G4cout << " Mom " << GammaDirection0 << G4endl;
166  G4cout << " Polarization " << GammaPolarization << G4endl;
167  G4cout << " MaterialPol. " << ElectronPolarization << G4endl;
168  G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
169  G4cout << " Log. Volume " << aLVolume->GetName() << G4endl;
170  G4cout << " Material " << aMaterial << G4endl;
171  }
172 
174  G4PhysicsVector * aVector=(*theAsymmetryTable)(midx);
175 
176  G4double asymmetry=0;
177  if (aVector) {
178  asymmetry = aVector->Value(GammaEnergy);
179  } else {
180  G4cout << " MaterialIndex " << midx << " is out of range \n";
181  asymmetry=0;
182  }
183 
184  // we have to determine angle between particle motion
185  // and target polarisation here
186  // circ pol * Vec(ElectronPol)*Vec(PhotonMomentum)
187  // both vectors in global reference frame
188 
189  G4double pol=ElectronPolarization*GammaDirection0;
190 
191  G4double polProduct = GammaPolarization.p3() * pol;
192  mfp *= 1. / ( 1. + polProduct * asymmetry );
193 
194  if (verboseLevel>=2) {
195  G4cout << " MeanFreePath: " << mfp / mm << " mm " << G4endl;
196  G4cout << " Asymmetry: " << asymmetry << G4endl;
197  G4cout << " PolProduct: " << polProduct << G4endl;
198  }
199  }
200 
201  return mfp;
202 }
203 
205  const G4Track& aTrack,
206  G4double previousStepSize,
208 {
209  // *** get unploarised mean free path from lambda table ***
210  G4double mfp = G4VEmProcess::PostStepGetPhysicalInteractionLength(aTrack, previousStepSize, condition);
211 
212 
214  // *** get asymmetry, if target is polarized ***
215  const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
216  const G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
217  const G4StokesVector GammaPolarization = aTrack.GetPolarization();
218  const G4ParticleMomentum GammaDirection0 = aDynamicGamma->GetMomentumDirection();
219 
220  G4Material* aMaterial = aTrack.GetMaterial();
221  G4VPhysicalVolume* aPVolume = aTrack.GetVolume();
222  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
223 
224  // G4Material* bMaterial = aLVolume->GetMaterial();
226 
227  const G4bool VolumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
228  G4StokesVector ElectronPolarization = polarizationManger->GetVolumePolarization(aLVolume);
229 
230  if (!VolumeIsPolarized || mfp == DBL_MAX) return mfp;
231 
232  if (verboseLevel>=2) {
233 
234  G4cout << " Mom " << GammaDirection0 << G4endl;
235  G4cout << " Polarization " << GammaPolarization << G4endl;
236  G4cout << " MaterialPol. " << ElectronPolarization << G4endl;
237  G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
238  G4cout << " Log. Volume " << aLVolume->GetName() << G4endl;
239  G4cout << " Material " << aMaterial << G4endl;
240  }
241 
243  G4PhysicsVector * aVector=(*theAsymmetryTable)(midx);
244 
245  G4double asymmetry=0;
246  if (aVector) {
247  asymmetry = aVector->Value(GammaEnergy);
248  } else {
249  G4cout << " MaterialIndex " << midx << " is out of range \n";
250  asymmetry=0;
251  }
252 
253  // we have to determine angle between particle motion
254  // and target polarisation here
255  // circ pol * Vec(ElectronPol)*Vec(PhotonMomentum)
256  // both vectors in global reference frame
257 
258  G4double pol=ElectronPolarization*GammaDirection0;
259 
260  G4double polProduct = GammaPolarization.p3() * pol;
261  mfp *= 1. / ( 1. + polProduct * asymmetry );
262 
263  if (verboseLevel>=2) {
264  G4cout << " MeanFreePath: " << mfp / mm << " mm " << G4endl;
265  G4cout << " Asymmetry: " << asymmetry << G4endl;
266  G4cout << " PolProduct: " << polProduct << G4endl;
267  }
268  }
269 
270  return mfp;
271 }
272 
273 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
274 
276 {
278 
280  G4bool isMaster = true;
281  const G4PolarizedCompton* masterProcess =
282  static_cast<const G4PolarizedCompton*>(GetMasterProcess());
283  if(masterProcess && masterProcess != this) { isMaster = false; }
284  if(isMaster) {
287  }
288  }
289 }
290 
291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
292 
293 
295 {
296  // *** build (unpolarized) cross section tables (Lambda)
298  if(buildAsymmetryTable && emModel) {
299  G4bool isMaster = true;
300  const G4PolarizedCompton* masterProcess =
301  static_cast<const G4PolarizedCompton*>(GetMasterProcess());
302  if(masterProcess && masterProcess != this) { isMaster = false; }
303  if(isMaster) { BuildAsymmetryTable(part); }
304  }
305 }
306 
307 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
308 
309 
311 {
312  // Access to materials
313  const G4ProductionCutsTable* theCoupleTable=
315  size_t numOfCouples = theCoupleTable->GetTableSize();
316  if(!theAsymmetryTable) { return; }
317  G4int nbins = LambdaBinning();
318  G4double emin = MinKinEnergy();
319  G4double emax = MaxKinEnergy();
320  G4PhysicsLogVector* aVector = 0;
321  G4PhysicsLogVector* bVector = 0;
322 
323  for(size_t i=0; i<numOfCouples; ++i) {
324  if (theAsymmetryTable->GetFlag(i)) {
325 
326  // create physics vector and fill it
327  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
328  // use same parameters as for lambda
329  if(!aVector) {
330  aVector = new G4PhysicsLogVector(emin, emax, nbins);
331  aVector->SetSpline(true);
332  bVector = aVector;
333  } else {
334  bVector = new G4PhysicsLogVector(*aVector);
335  }
336 
337  for (G4int j = 0; j <= nbins; ++j ) {
338  G4double energy = bVector->Energy(j);
339  G4double tasm=0.;
340  G4double asym = ComputeAsymmetry(energy, couple, part, 0., tasm);
341  bVector->PutValue(j,asym);
342  }
343 
345  }
346  }
347 }
348 
349 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
350 
351 
353  const G4MaterialCutsCouple* couple,
354  const G4ParticleDefinition& aParticle,
355  G4double cut,
356  G4double & tAsymmetry)
357 {
358  G4double lAsymmetry = 0.0;
359  tAsymmetry=0;
360 
361  //
362  // calculate polarized cross section
363  //
364  G4ThreeVector thePolarization=G4ThreeVector(0.,0.,1.);
365  emModel->SetTargetPolarization(thePolarization);
366  emModel->SetBeamPolarization(thePolarization);
367  G4double sigma2=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
368 
369  //
370  // calculate unpolarized cross section
371  //
372  thePolarization=G4ThreeVector();
373  emModel->SetTargetPolarization(thePolarization);
374  emModel->SetBeamPolarization(thePolarization);
375  G4double sigma0=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
376 
377  // determine assymmetries
378  if (sigma0 > 0.) {
379  lAsymmetry = sigma2/sigma0-1.;
380  }
381  return lAsymmetry;
382 }
383 
384 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4double condition(const G4ErrorSymMatrix &m)
const G4VProcess * GetMasterProcess() const
Definition: G4VProcess.hh:538
const G4ThreeVector & GetPolarization() const
static const double MeV
Definition: G4SIunits.hh:193
static void SetPhysicsVector(G4PhysicsTable *physTable, size_t idx, G4PhysicsVector *vec)
G4double MaxKinEnergy() 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)
G4VEmModel * EmModel(G4int index=1) const
static G4PolarizationManager * GetInstance()
virtual void PreparePhysicsTable(const G4ParticleDefinition &)
void SetSplineFlag(G4bool val)
G4double MaxKinEnergy() const
G4PolarizedComptonModel * emModel
void BuildPhysicsTable(const G4ParticleDefinition &)
G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
void SetStartFromNullFlag(G4bool val)
int G4int
Definition: G4Types.hh:78
static G4PhysicsTable * PreparePhysicsTable(G4PhysicsTable *physTable)
G4int LambdaBinning() const
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
void SetSpline(G4bool)
G4double p3() const
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:707
void SetMinKinEnergyPrim(G4double e)
virtual G4double GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)
void SetEmModel(G4VEmModel *, G4int index=1)
size_t CurrentMaterialCutsCoupleIndex() const
G4GLOB_DLL std::ostream G4cout
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:483
void PutValue(size_t index, G4double theValue)
G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
G4double MinKinEnergy() const
void PreparePhysicsTable(const G4ParticleDefinition &)
G4double Energy(size_t index) const
G4double Value(G4double theEnergy, size_t &lastidx) const
G4Material * GetMaterial() const
static G4ProductionCutsTable * GetProductionCutsTable()
G4PolarizedCompton(const G4String &processName="pol-compt", G4ProcessType type=fElectromagnetic)
G4double ComputeAsymmetry(G4double energy, const G4MaterialCutsCouple *couple, const G4ParticleDefinition &particle, G4double cut, G4double &tAsymmetry)
G4LogicalVolume * GetLogicalVolume() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
void SetBeamPolarization(const G4ThreeVector &pBeam)
void AddEmModel(G4int, G4VEmModel *, const G4Region *region=0)
G4double energy(const ThreeVector &p, const G4double m)
void SetSecondaryParticle(const G4ParticleDefinition *p)
static G4EmParameters * Instance()
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:777
G4bool GetFlag(size_t i) const
double G4double
Definition: G4Types.hh:76
static G4PhysicsTable * theAsymmetryTable
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:714
G4ForceCondition
G4ThreeVector G4ParticleMomentum
G4double MinKinEnergy() const
#define DBL_MAX
Definition: templates.hh:83
static const double mm
Definition: G4SIunits.hh:102
const G4ThreeVector & GetVolumePolarization(G4LogicalVolume *lVol) const
void SetTargetPolarization(const G4ThreeVector &pTarget)
virtual void PrintInfo()
void SetModel(const G4String &name)
G4ProcessType